python 数据上的蒙特卡罗拟合

如何解决python 数据上的蒙特卡罗拟合

我编写了一个蒙特卡罗模拟来拟合 49 个具有非对称误差条的数据点。由于误差在两个轴上都是不对称的，我不能简单地使用 scipy.optimize.curve_fit 模块。这是我的基本方法：

在某个数据点使用具有最大概率的三角概率分布分布，从置信水平（误差范围）内生成一个包含 1000 个随机数的列表。现在我有一个维度列表 [49*1000]。
将此列表从 [491000] 转换为 [100049]。我这样做是为了得到一个包含 1000 个样本的数据集，这些样本包含 49 个在误差范围内的点。
使用 scipy 的 curve_fit 分别拟合这 1000 个样本，并在函数 y=m*x+c 中找到自由参数（c 是自由参数，我已经知道斜率 m。）
使用 sklearn.metrics.mean_squared_error 模块在这 1000 个样本中的每一个中查找均方误差。
找到具有最小均方值的索引，并使用该索引处的 popt 值（c 参数）绘制拟合图。

工作代码：这是我的代码：

trials=int(1e4)
    
#xdata: most probable x value
#ydata: most probable y value
#xerr6low: lower bound on x error
#xerr6up: upper bound on x error
#yerr6low: lower bound on y error
#yerr6up: upper bound on y error

#generating random number using triangular distribution weighted with highest probability at xdata/ydata:
xarray1=[0]*len(obs6xr)
yarray1=[0]*len(obs6xr)
for i in range(0,len(obs6xr)):
    xarray1[i]=np.random.triangular(xdata[i]-xerr6low[i],xdata[i],xdata[i]+xerr6up[i],trials)
    yarray1[i]=np.random.triangular(ydata[i]-yerr6low[i],ydata[i],ydata[i]+yerr6up[i],trials)

xarray=[list(x) for x  in zip(*xarray1)]
yarray=[list(x) for x  in zip(*yarray1)]


def func(x,c):
    #return (np.log10(a)+(b*(x-12))+np.log10(10**(8+c)))
    #return ((x/1e12)**a)*(b)*(10**(8+c))
    m = 1.65
    mx = [element * m for element in x]
    y = [j+c for j in mx]
    return y

#Fit for the parameters a,b,c of the function func:
popt=np.zeros(trials)
pcov=np.zeros(trials)
print(len(xarray),len(yarray),len(popt))
for i in tqdm(range (trials),desc='Optimizing'):
    popt[i],pcov[i] = curve_fit(func,xarray[i],np.array(yarray[i]))

MSE=np.zeros(trials)
for i in tqdm(range (trials),desc='Calculating MSE'):
    MSE[i]=mean_squared_error(np.array(yarray[i]),func(np.array(xarray[i]),popt[i]))

minimum_MSE=np.amin(MSE)
index_MSE=np.where(MSE == np.amin(MSE))
print('minimum MSE = ',minimum_MSE,'at index = ',index_MSE)
print(popt[index_MSE])

def MbhthShimasakupropc(x,c):
    Mbhth = (10**(c))*(x**1.65)
    return Mbhth

plt.figure(figsize=(5,5))
x=np.logspace(11,14,trials)
plt.loglog(obs6x,MbhthShimasakuprop(obs6x,popt[index_MSE]))
plt.xscale('log')
plt.yscale('log')
m=np.logspace(np.log10(4e10),np.log10(3e14),1000)
plt.errorbar(obs0x,obs0y,xerr=asymmetric_errorx0,yerr=asymmetric_errory0,fmt='o',color='black',markersize='2.5',ecolor='black',capsize=2,elinewidth=1)
plt.errorbar(obs6x,obs6y,xerr=asymmetric_errorx6,yerr=asymmetric_errory6,color='red',ecolor='red',elinewidth=1)
plt.loglog(m,MbhthShimasaku(m,0),linestyle=':',label='Shimasaku-Ferrarese z=0')

plt.xlim(4e10,3e14)
plt.ylim(1e6,1e11)
plt.legend(['Monte-Carlo Fitting','local relation','z=0','z~6'])
plt.show()

结果：

正如我们所见，这段代码运行良好。但是如果我将函数更改为 2 个参数 (a&b)， func(x,a,b) 并应用相同的代码并进行一些小的调整，代码失败得很惨。

不工作的代码：

trials=int(1e3)

#generating random number:
xarray1=[0]*len(obs6xr)
yarray1=[0]*len(obs6xr)
for i in range(0,trials)

xarray=[list(x) for x  in zip(*xarray1)]
yarray=[list(x) for x  in zip(*yarray1)]

def func(x,b):
    y=a*x+b
    return y    

#Fit for the parameters a,c of the function func:
popt=[0]*(trials)
pcov=[0]*(trials)
print(len(xarray),len(popt))
for i in tqdm(range (trials)):
    popt[i],np.array(xarray[i]),np.array(yarray[i]))

   
MSE=np.zeros(trials)
for i in tqdm(range (trials)):
    MSE[i]=mean_squared_error(np.array(yarray[i]),np.array(popt[i])[0],np.array(popt[i])[1]))

minimum_MSE=np.amin(MSE)
index_MSE=np.where(MSE == np.amin(MSE))
print('minimum MSE = ',index_MSE)


def MbhthShimasakupropmc(x,m,c):
    Mbhth = (10**(c))*(x**m)
    return Mbhth

plt.figure(figsize=(5,MbhthShimasakupropmc(obs6x,*popt[[index_MSE][0][0][0]]))
plt.xscale('log')
plt.yscale('log')
m=np.logspace(np.log10(4e10),label='Shimasaku-Ferrarese z=0')

plt.legend(['Monte-Carlo Fitting','z~6'])
plt.show()

结果：

我不知道我做错了什么。有人可以帮我调试问题吗。

python 数据上的蒙特卡罗拟合

如何解决python 数据上的蒙特卡罗拟合

相关推荐