与 python curve_fit 的灾难性契合
catastrophic fit with python curve_fit
我需要用指数模型拟合数据(x 轴:sigma,y 轴:Mbh)。这是我的代码:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
#define my data
Mbh = np.array([1.8e6,2.5e6,4.5e7,3.7e7,4.4e7,1.5e7,1.4e7,4.1e7, 1.0e9,2.1e8,1.0e8,1.0e8,1.6e7,1.9e8,3.9e7,5.2e8,3.1e8,3.0e8,7.0e7,1.1e8,3.0e9,5.6e7,7.8e7,2.0e9,1.7e8,1.4e7,2.4e8,5.3e8,3.3e8,3.5e6,2.5e9])
sigma = np.array([103,75,160,209,205,151,175,140,230,205,145,206,143,182,130,315,242,225,186,190,375,162,152,385,177,90,234,290,266,67,340])
#define my model to fit
def Mbh02(alpha, sigma, beta):
return alpha * np.exp(beta*sigma);
#calculate the fit parameter:
#for second model
popt02, pcov02 = curve_fit(Mbh02, sigma, Mbh, p0=[1, 0.058])
print(f'Parameter of the second function : {popt02}')
sigma_plot = [103,75,160,209,205,151,175,140,230,205,145,206,143,182,130,315,242,225,186,190,375,162,152,385,177,90,234,290,266,67,340]
sigma_plot.sort()
sigma_plot = np.array(sigma_plot)
#plot model with data with
plt.figure(figsize=(6,6))
plt.scatter(sigma, Mbh * 1e-9, marker = '+', color ='black', label = 'Data')
plt.plot(sigma_plot , Mbh02(alpha = popt02[0], sigma = sigma_plot, beta = popt02[1]) * 1e-9, color='orange', ls ='-', label ='2. fit')
plt.ylabel(r'$M_{BH}$ in $M_\odot *10^9$ unit', fontsize=16)
plt.xlabel(r'$\sigma$', fontsize=16)
# plt.ylim(-1,10)
plt.title('Plot of the black hole mass $M_{BH}$ \nagainst the velocity dispersion $\sigma$ \nfor different elliptical galaxies', fontsize=18)
plt.grid(True)
plt.legend()
plt.show()
我得到以下参数
:
print(popt01) = [16.13278858 0.91788691]
看起来:
如果我尝试手动查找参数,并绘制它们:
plt.plot(sigma_plot , (1 * np.exp(0.058 * sigma_plot)) * 1e-9, ls ='--', label ='2. fit manual')
我得到了下面的情节,它好多了:
有什么问题吗?为什么 curve_fit 不工作并给出这样的参数?
Assumes ydata = f(xdata, *params) + eps
因此,如果您更改函数定义,使 x 数据在您的函数中位于第一个,它将起作用:
def Mbh02(sigma, alpha, beta):
return alpha * np.exp(beta*sigma);
# Rest of code
plt.plot(sigma_plot , Mbh02(sigma_plot, *popt02) * 1e-9, color='orange', ls ='-')
您是否尝试过使用线性拟合而不是拟合 exp 来拟合 log(Mbh)。直接建模?这通常会带来很大的稳定性。
import numpy as np
import matplotlib.pyplot as plt
Mbh = np.array([1.8e6,2.5e6,4.5e7,3.7e7,4.4e7,1.5e7,1.4e7,4.1e7, 1.0e9,2.1e8,1.0e8,1.0e8,1.6e7,1.9e8,3.9e7,5.2e8,3.1e8,3.0e8,7.0e7,1.1e8,3.0e9,5.6e7,7.8e7,2.0e9,1.7e8,1.4e7,2.4e8,5.3e8,3.3e8,3.5e6,2.5e9])
sigma = np.array([103,75,160,209,205,151,175,140,230,205,145,206,143,182,130,315,242,225,186,190,375,162,152,385,177,90,234,290,266,67,340])
plt.figure(2)
plt.plot(sigma,Mbh,'.')
lnMbh= np.log(Mbh)
p = np.polyfit(sigma,lnMbh,1)
plt.plot(sigma, np.exp(np.polyval(p,sigma)),'*')
alpha = np.log(p[0])
beta = p[1]
我需要用指数模型拟合数据(x 轴:sigma,y 轴:Mbh)。这是我的代码:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
#define my data
Mbh = np.array([1.8e6,2.5e6,4.5e7,3.7e7,4.4e7,1.5e7,1.4e7,4.1e7, 1.0e9,2.1e8,1.0e8,1.0e8,1.6e7,1.9e8,3.9e7,5.2e8,3.1e8,3.0e8,7.0e7,1.1e8,3.0e9,5.6e7,7.8e7,2.0e9,1.7e8,1.4e7,2.4e8,5.3e8,3.3e8,3.5e6,2.5e9])
sigma = np.array([103,75,160,209,205,151,175,140,230,205,145,206,143,182,130,315,242,225,186,190,375,162,152,385,177,90,234,290,266,67,340])
#define my model to fit
def Mbh02(alpha, sigma, beta):
return alpha * np.exp(beta*sigma);
#calculate the fit parameter:
#for second model
popt02, pcov02 = curve_fit(Mbh02, sigma, Mbh, p0=[1, 0.058])
print(f'Parameter of the second function : {popt02}')
sigma_plot = [103,75,160,209,205,151,175,140,230,205,145,206,143,182,130,315,242,225,186,190,375,162,152,385,177,90,234,290,266,67,340]
sigma_plot.sort()
sigma_plot = np.array(sigma_plot)
#plot model with data with
plt.figure(figsize=(6,6))
plt.scatter(sigma, Mbh * 1e-9, marker = '+', color ='black', label = 'Data')
plt.plot(sigma_plot , Mbh02(alpha = popt02[0], sigma = sigma_plot, beta = popt02[1]) * 1e-9, color='orange', ls ='-', label ='2. fit')
plt.ylabel(r'$M_{BH}$ in $M_\odot *10^9$ unit', fontsize=16)
plt.xlabel(r'$\sigma$', fontsize=16)
# plt.ylim(-1,10)
plt.title('Plot of the black hole mass $M_{BH}$ \nagainst the velocity dispersion $\sigma$ \nfor different elliptical galaxies', fontsize=18)
plt.grid(True)
plt.legend()
plt.show()
我得到以下参数 :
print(popt01) = [16.13278858 0.91788691]
看起来:
如果我尝试手动查找参数,并绘制它们:
plt.plot(sigma_plot , (1 * np.exp(0.058 * sigma_plot)) * 1e-9, ls ='--', label ='2. fit manual')
我得到了下面的情节,它好多了:
有什么问题吗?为什么 curve_fit 不工作并给出这样的参数?
Assumes ydata = f(xdata, *params) + eps
因此,如果您更改函数定义,使 x 数据在您的函数中位于第一个,它将起作用:
def Mbh02(sigma, alpha, beta):
return alpha * np.exp(beta*sigma);
# Rest of code
plt.plot(sigma_plot , Mbh02(sigma_plot, *popt02) * 1e-9, color='orange', ls ='-')
您是否尝试过使用线性拟合而不是拟合 exp 来拟合 log(Mbh)。直接建模?这通常会带来很大的稳定性。
import numpy as np
import matplotlib.pyplot as plt
Mbh = np.array([1.8e6,2.5e6,4.5e7,3.7e7,4.4e7,1.5e7,1.4e7,4.1e7, 1.0e9,2.1e8,1.0e8,1.0e8,1.6e7,1.9e8,3.9e7,5.2e8,3.1e8,3.0e8,7.0e7,1.1e8,3.0e9,5.6e7,7.8e7,2.0e9,1.7e8,1.4e7,2.4e8,5.3e8,3.3e8,3.5e6,2.5e9])
sigma = np.array([103,75,160,209,205,151,175,140,230,205,145,206,143,182,130,315,242,225,186,190,375,162,152,385,177,90,234,290,266,67,340])
plt.figure(2)
plt.plot(sigma,Mbh,'.')
lnMbh= np.log(Mbh)
p = np.polyfit(sigma,lnMbh,1)
plt.plot(sigma, np.exp(np.polyval(p,sigma)),'*')
alpha = np.log(p[0])
beta = p[1]