здесь нужна ваша помощь!
Я генерирую зашумленные данные Лоренца и подгоняю их с помощью функции Лоренца, используя scipy.optimize.curve_fit . Моя цель состоит в том, чтобы создать посадки с таким же уровнем шума (равномерный шум, как показано в коде ниже), но с увеличением количества точек данных, которые установлены. Мне удается это сделать, и результаты показаны ниже (где количество точек данных, которые установлены go от 10 до 100 с шагом 10).
Теперь я хочу выполнить подгонку для каждого количества точек данных 1000 раз (на этот раз с 10 до 100 точек с шагом 1); так что я могу оценить среднюю ошибку (определяемую как разницу между истинным и оценочным значением по модели) по мере того, как меняется число подгоненных точек. Когда я отображаю предполагаемую ошибку как функцию точек данных, она, кажется, колеблется вверх и вниз (а не плавно спускаясь по мере увеличения числа точек). Любая идея, почему это? Любая помощь будет принята с благодарностью. Код прилагается, а изображения приведены после кода в виде ссылок:
#%%estimation of parameters and 95 % confidence bounds
# Initial Lorentzian Parameters
y_0 = 1 #half width at half maximum (HWHM)
x_0 = 0 #peak of distribution
I= 1 #amplitude
#Lorentzian Model
def Lorentzian(x, x_0, y_0, I):
return (I*y_0**2)/((x_array-x_0)**2+y_0**2)
#declare arrays to be filled later
peak_center = []
peak_center_err = []
hwhm = []
hwhm_err = []
amp = []
amp_err = []
j=0
while j <1000: #running simulation 1000 times
j=j+1
for i in range(10,101,1): #from 10 to 100 in steps of 1 (testing 90 point conditions)
#linearly spaced x arrays from -10 to +10 (containing from 10 to 100 points in steps of 1)
x_array = np.linspace(-10,10,i)
#defining lorentzian
y_array_lorentzian = ((I*y_0**2)/((x_array-x_0)**2+y_0**2))/I #normalised
# creating some noise to add to the the y-axis data
#Uniform noise (positive)
#Return random floats in the half-open interval [0.0, 1.0) from continous unifrom distribution
noise = np.random.ranf(i)*0.01
y_array_lorentzian += noise
#fitting
#Use non-linear least squares to fit a function, Lorentzian, to data.
popt_lorentz, pcov_lorentz = scipy.optimize.curve_fit(Lorentzian, x_array, y_array_lorentzian,p0=[1,1,1],
absolute_sigma=False)
perr_lorentz = np.sqrt(np.diag(pcov_lorentz))
#storing optimal values and error
peak_center.append(popt_lorentz[0])
peak_center_err.append(abs(popt_lorentz[0]-x_0))
hwhm.append(popt_lorentz[1])
hwhm_err.append(abs(popt_lorentz[1]-y_0))
amp.append(popt_lorentz[2])
amp_err.append(abs(popt_lorentz[2]-I))
#reshaping arrays
#error of peak center
peak_center_err=np.array(peak_center_err)
peak_center_err = peak_center_err.reshape(1000,91)
#error of hwhm
hwhm_err = np.array(hwhm_err)
hwhm_err = hwhm_err.reshape(1000,91)
#error of amp
amp_err = np.array(amp_err)
amp_err = amp_err.reshape(1000,91)
##note: 1000 rows --> number of simulations (j) AND 90 columns --> increasing number of points (i)
#this means that first column contains all estimations with 10 data points, second column all estimations with 11 points etc
#find mean of error and sd of error for parameters
mean_err_peak_center =[]
sd_err_peak_center=[]
mean_err_hwhm=[]
sd_err_hwhm=[]
mean_err_amp = []
sd_err_amp = []
for k in range(91):
#mean
mean_err_peak_center.append(peak_center_err[:,k].mean(axis=0))
mean_err_hwhm.append(hwhm_err[:,k].mean(axis=0))
mean_err_amp.append(amp_err[:,k].mean(axis=0))
#standard deviation
sd_err_peak_center.append(peak_center_err[:,k].std())
sd_err_hwhm.append(hwhm_err[:,k].std())
sd_err_amp.append(amp_err[:,k].std())
#95 % bounds of error
ci_peak_center = 2*np.array(sd_err_peak_center)
ci_hwhm = 2*np.array(sd_err_hwhm)
ci_amp = 2*np.array(sd_err_amp)
#plot of mean error estimation and 95 % bounds of error
points = np.linspace(10,100,91)
fig=plt.figure(figsize=(12,4))
ax=fig.add_subplot(1,3,1)
ax.plot(points,mean_err_peak_center, c= 'k', marker = 'o', ms=3)
ax.fill_between(points, (mean_err_peak_center-ci_peak_center), (mean_err_peak_center+ci_peak_center), color='r', alpha=.1)
ax.set_title('Mean Error $x_0$ (Peak center)')
ax.set_xlabel('Number of Points')
ax.set_ylabel('$|x_{0 est} - x_{0 true}|$ [GHz]')
ax.legend(['Mean Error', '95 % Bounds'], frameon=False)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
#plt.xticks(ticks=[i for i in range(10,110,5)], labels=[i for i in range(10,110,5)])
ax1= fig.add_subplot(1,3,2)
ax1.plot(points,mean_err_hwhm, linewidth = 0.5, c= 'k', ls = '-', marker = 'o', ms=3)
ax1.fill_between(points, (mean_err_hwhm-ci_hwhm), (mean_err_hwhm+ci_hwhm), color='r', alpha=.1)
ax1.set_title('Mean Error $y_0$ (HWHM)')
ax1.set_xlabel('Number of Points')
ax1.set_ylabel('$|y_{0 est} - y_{0 true}|$ [GHz]')
ax1.legend(['Mean Error', '95 % Bounds'], frameon=False)
ax1.spines['right'].set_visible(False)
ax1.spines['top'].set_visible(False)
#plt.xticks(ticks=[i for i in range(10,110,5)], labels=[i for i in range(10,110,5)])
ax2= fig.add_subplot(1,3,3)
ax2.plot(points,mean_err_amp, linewidth = 0.5, c= 'k', ls = '-', marker = 'o', ms=3)
ax2.fill_between(points, (mean_err_amp-ci_amp), (mean_err_amp+ci_amp), color='r', alpha=.1)
ax2.set_title('Mean Error $I$ (Amplitude)')
ax2.set_xlabel('Number of Points')
ax2.set_ylabel('$|I_{est} - I_{true}|$ [a.u.]')
ax2.legend(['Mean Error', '95 % Bounds'], frameon=False)
ax2.spines['right'].set_visible(False)
ax2.spines['top'].set_visible(False)
#plt.xticks(ticks=[i for i in range(10,110,5)], labels=[i for i in range(10,110,5)])
fig.subplots_adjust(hspace=0.5, wspace = 0.5)
fig.tight_layout()
соответствует лоренцу
расчетная ошибка 1000 симуляций