拟合函数参数和误差--最小均方根

拟合函数参数和误差——方法（1）leastsq

python的Scipy.optimize中有函数可以用来拟合函数，可以用来求参数和误差。需要有一组数据X，假设其为高斯分布，求其分布参数$\mu ,\sigma$。可以先有scipy的leastsq计算其初步的参数，然后再由lmfit模块求更精确的参数以及其不确定度。代码为：

import math
import numpy as np
import lmfit
import matplotlib.pyplot as plt
from scipy.optimize import leastsq



def Gauss(t,para):
    A,sig,tp = para
    #print(t)
    return A*math.e**(-((t-tp)**2)/(2*(sig**2)))
 
 
###需要拟合的函数func及误差error###
def ReGauss(para,pul,t):
   ### 误差
    #print("ERR:",(pul - Gauss(t,para))**2)
    return (pul - Gauss(t,para))
def residual(p):
    v = p.valuesdict()
    return v['A']*math.e**(-((x-v['mean'])**2)/(2*(v['sig']**2))) -y
    

bwith=1.5

np.random.seed(1000)
alpha = np.random.normal(-2,3,100) #存成二进制序列、方便读取

num_bin = 24


n_hist,edges = np.histogram(alpha,bins=num_bin,range=(-8,4))
num_hist = np.array(n_hist)

alpha_m = np.mean(alpha)
alpha_std = np.std(alpha)

alpha_x = np.arange(-8,4,0.5)

p_init = [1,alpha_std,alpha_m]
para,pcov,infodict,errmsg,success = leastsq(ReGauss,p_init,args=(n_hist,alpha_x),full_output=1)


## 计算 拟合参数的不确定度
para_lmft =  lmfit.Parameters()
para_lmft.add_many(('A',para[0]),('sig',para[1]),('mean',para[2]))
x =  np.arange(-8,4,0.5)
y = n_hist
mi = lmfit.minimize(residual, para_lmft,method='nelder', nan_policy='omit')
lmfit.printfuncs.report_fit(mi.params, min_correl=0.5)


res = lmfit.minimize(residual, method='emcee', nan_policy='omit', burn=300, steps=1000,thin=20,
                     params=mi.params, is_weighted=False, progress=False)

highest_prob = np.argmax(res.lnprob)
hp_loc = np.unravel_index(highest_prob, res.lnprob.shape)
mle_soln = res.chain[hp_loc]
for i, par in enumerate(para_lmft):
    para_lmft[par].value = mle_soln[i]


print('\nMaximum Likelihood Estimation from emcee       ')
print('-------------------------------------------------')
print('Parameter  MLE Value   Median Value   Uncertainty')
fmt = '  {:5s}  {:11.5f} {:11.5f}   {:11.5f}'.format

#print(para_lmft.items(name))

para_last = []
err = []
for name, param in para_lmft.items():
    para_last.append(param.value)
    err.append(res.params[name].stderr)
    print(fmt(name, param.value, res.params[name].value,res.params[name].stderr))

print("参数A,sig,mean为：",para_last)
print("参数A,sig,mean不确定度为:",err)


alpha_x = np.linspace(-12,8,10000)
alpha_y = Gauss(alpha_x,para_last)

print("Mean alpha",alpha_m)
plt.figure(figsize=(12,9),dpi=100)

ax = plt.subplot()

plt.tick_params(labelsize=13.5)
ax.tick_params(which='both',top=True,bottom=True,left=True,right=True)
ax.spines['bottom'].set_linewidth(bwith)
ax.spines['left'].set_linewidth(bwith)
ax.spines['top'].set_linewidth(bwith)
ax.spines['right'].set_linewidth(bwith)
#axes.xaxis.set_minor_locator( MultipleLocator(0.1) )


plt.tick_params(which='both',top=True,bottom='on',left='on',right='on',labelsize=12)
plt.hist(alpha,bins=24,color="w",edgecolor="black",lw=2,alpha=0.8,range=(-8,4))
plt.plot(alpha_x,alpha_y,color="r",lw=2)

plt.ylabel("Number of sample",fontsize=15)
plt.xlabel("Some Para X",fontsize=15)
plt.show()
  

结果图为：

---

End