使用最佳拟合生成 chi_squared
Generating chi_squared using best fit
我正在生成一些假数据集的 1000 次迭代,并使用 curve_fit 找到最佳拟合模型(具有一定均值、偏移量、amp...的高斯模型)而不是猜测我的平均值拟合模型,使用 curve_fit 以便我可以获得更好的 chi2 值。
但是现在输出很奇怪,甚至数据和模型也不一致。
此外,如果我绘制 chi2 值的直方图,
true_mean = 5 #Assume that we know the perfect fitting model is gaussian
#with mean value 5. And use try_mean with different mean values to see how
#does the chi2 behave
true_amp = 100
true_wid = 2
true_offset =10
x_values = np.array([i for i in np.arange(0,10,0.4)])
exact_y_values = np.array([true_offset + true_amp*
np.exp(-((i-true_mean)/true_wid)**2)
for i in x_values])
try_mean = 4.8 # Notice the data is generated centered at 5;
# by comparing it to a 4.8 we expect disagreement.
try_amp = 100
try_wid = 2
try_offset =10
try_these_y_values = np.array([try_offset + try_amp*
np.exp(-((i-try_mean)/try_wid)**2)
for i in x_values])
def func (x_values,offset,amp,mean,wid):
return (offset + amp*np.exp(-((i-mean)/wid)**2))
#Excercise 2
#def func (x_values,offset,amp,mean,wid): return (offset + amp*np.exp(-((i-mean)/wid)**2))
chi2_fit=np.zeros(1000)
for i in range (1000):
fake_exp_y_values = np.array([np.random.poisson(y)
for y in exact_y_values])
p01=[true_offset,true_amp,true_mean,true_wid]
[popt,pcov]=opt.curve_fit(func,x_values, fake_exp_y_values,p01)
y_values_fit=np.array([popt[0]+ popt[1]
*np.exp(
-((x_values-popt[2])/popt[3])**2)])
residuals=fake_exp_y_values-y_values_fit
y_err=np.clip(np.sqrt(fake_exp_y_values),1,9999)
pulls=residuals/y_err
chi2_fit[i]=np.sum(pulls**2)
plt.hist(chi2_fit)
plt.scatter(x_values,exact_y_values,color='k',ls='--',
label='true model')
plt.scatter(x_values,y_values_fit,color='r')
plt.errorbar(x_values,y_values_fit,yerr=y_err)
为了剧情合理应该修改什么
问题是你函数的定义和用法 i
:
def func (x_values,offset,amp,mean,wid):
return (offset + amp*np.exp(-((i-mean)/wid)**2))
for i in range (1000):
...
[popt,pcov]=opt.curve_fit(func,x_values, fake_exp_y_values,p01)
...
更正版本:
def func (x_values,offset,amp,mean,wid):
return (offset + amp*np.exp(-((x_values-mean)/wid)**2))
下次还请postMRE。这个不是最小的,也不是可复制的(我必须在最后添加 y_values_fit = y_values_fit.flatten()
才能使其工作)。
我正在生成一些假数据集的 1000 次迭代,并使用 curve_fit 找到最佳拟合模型(具有一定均值、偏移量、amp...的高斯模型)而不是猜测我的平均值拟合模型,使用 curve_fit 以便我可以获得更好的 chi2 值。
但是现在输出很奇怪,甚至数据和模型也不一致。
此外,如果我绘制 chi2 值的直方图,
true_mean = 5 #Assume that we know the perfect fitting model is gaussian
#with mean value 5. And use try_mean with different mean values to see how
#does the chi2 behave
true_amp = 100
true_wid = 2
true_offset =10
x_values = np.array([i for i in np.arange(0,10,0.4)])
exact_y_values = np.array([true_offset + true_amp*
np.exp(-((i-true_mean)/true_wid)**2)
for i in x_values])
try_mean = 4.8 # Notice the data is generated centered at 5;
# by comparing it to a 4.8 we expect disagreement.
try_amp = 100
try_wid = 2
try_offset =10
try_these_y_values = np.array([try_offset + try_amp*
np.exp(-((i-try_mean)/try_wid)**2)
for i in x_values])
def func (x_values,offset,amp,mean,wid):
return (offset + amp*np.exp(-((i-mean)/wid)**2))
#Excercise 2
#def func (x_values,offset,amp,mean,wid): return (offset + amp*np.exp(-((i-mean)/wid)**2))
chi2_fit=np.zeros(1000)
for i in range (1000):
fake_exp_y_values = np.array([np.random.poisson(y)
for y in exact_y_values])
p01=[true_offset,true_amp,true_mean,true_wid]
[popt,pcov]=opt.curve_fit(func,x_values, fake_exp_y_values,p01)
y_values_fit=np.array([popt[0]+ popt[1]
*np.exp(
-((x_values-popt[2])/popt[3])**2)])
residuals=fake_exp_y_values-y_values_fit
y_err=np.clip(np.sqrt(fake_exp_y_values),1,9999)
pulls=residuals/y_err
chi2_fit[i]=np.sum(pulls**2)
plt.hist(chi2_fit)
plt.scatter(x_values,exact_y_values,color='k',ls='--',
label='true model')
plt.scatter(x_values,y_values_fit,color='r')
plt.errorbar(x_values,y_values_fit,yerr=y_err)
为了剧情合理应该修改什么
问题是你函数的定义和用法 i
:
def func (x_values,offset,amp,mean,wid):
return (offset + amp*np.exp(-((i-mean)/wid)**2))
for i in range (1000):
...
[popt,pcov]=opt.curve_fit(func,x_values, fake_exp_y_values,p01)
...
更正版本:
def func (x_values,offset,amp,mean,wid):
return (offset + amp*np.exp(-((x_values-mean)/wid)**2))
下次还请postMRE。这个不是最小的,也不是可复制的(我必须在最后添加 y_values_fit = y_values_fit.flatten()
才能使其工作)。