计算科学 - 为什么我的 curve_fit 没有生成协方差矩阵和未知变量的正确值？ - 吾爱随笔录

为什么我的 curve_fit 没有生成协方差矩阵和未知变量的正确值？

计算科学 Python scipy 曲线拟合

2021-12-18 11:00:02

我正在尝试将超新星数据拟合到一个scipy.curve_fit函数中。但是，当我的代码运行时，给出的未知变量的值popt是准确的。也没有产生协方差矩阵。该代码没有为未知变量提供正确的值 $K_1$ , $K_2$ , $\alpha$ 和 $\beta$ . 我使用的适合我的模型如下：

f = K_{1} ((1.39 / 5)^{α}) (t^{β}) (e^{- (K_{2} (1.39 / 5)^{- 2.1} t^{- 3}})

$f = K_1((1.39/5)^\alpha) (t^\beta) (e^{-(K_2(1.39/5)^{-2.1} t^{-3}})\,$

下面是我的数据集，其中年月日之后的第二列作为 t，第四列作为通量密度，第五列（最后一列）作为 yerr。

有人可以告诉如何在这段代码中产生协方差矩阵吗？另外，我的代码最初的猜测是什么？

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


t=np.genfromtxt("data_1390.dat",dtype=np.float64,usecols=3)
y=np.genfromtxt("data_1390.dat",dtype=np.float64,usecols=5)
yerr=np.genfromtxt("data_1390.dat",dtype=np.float64 ,usecols=6) 


def FFA(t,K1,K2,alpha,beta):
        CSM_unif=K2*((1.39/5)**(-2.1))*(t/1)**(-3)
        
        return K1*((1.39/5)**(alpha))*((t/1)**beta)*np.exp(-CSM_unif)
        
fig,ax1=plt.subplots()

plt.errorbar(t,y,yerr=yerr,color='blue',fmt='o-',elinewidth=3,capsize=2,barsabove=True)
#sigma=[132899.651,-261.945481,-261.945481,14.3752141]
popt,pcov=curve_fit(FFA,t,y,p0=[1.02,5,1.151,1.151], sigma=yerr, maxfev=2000)


K1,K2,alpha,beta=popt
print("value of K1 = ",K1)
print("value of K2 = ",K2)
print("value of beta = ",beta)
print("value of alpha = ",alpha)

ax1.set_xscale("log", nonpositive='clip')
ax1.set_yscale("log", nonpositive='clip')

ax1.plot(t,y,FFA(t,*popt),color='r')

ax1.set_xlabel('time (in s)')
ax1.set_ylabel('spectral flux density (in Jansky)')
ax1.set(title='Light Curve')

ax1.set_ylim(bottom=0.1)
ax1.set_xlim(left=10)
fig.tight_layout()
plt.show()

以下是我的警告信息和输出；

Warning (from warnings module):
  File "C:\Users\HP\AppData\Local\Programs\Python\Python39\lib\site-packages\scipy\optimize\minpack.py", line 833
    warnings.warn('Covariance of the parameters could not be estimated',
OptimizeWarning: Covariance of the parameters could not be estimated
value of K1 =  579.2356457671407
value of K2 =  243.98216449823164
value of beta =  1.1217015399194465
value of alpha =  444.59686079559594

我的最佳拟合曲线。曲线不正确，因为弯曲应该更高。最大值应该更高。

1个回答

问题似乎是缩放问题之一。当我添加函数的 jacobian 时，出现了溢出警告。因此，我将数据除以它们的最大值并且它起作用了。以下是代码。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

t = np.array([33.90, 76.95, 166.65, 302.15, 330.11, 429.82, 533.59, 638.19])
t = t/t.max()
y = np.array([0.25, 1.81, 8.32, 11.60, 12.18, 10.12, 9.44, 5.81])
y = y/y.max()


def FFA(t, K1, K2, beta, delta):
        CSM_unif = K2 * t**delta
        return K1 * t**beta * np.exp(-CSM_unif)
    
    
def grad(t, K1, K2, beta, delta):
    g1 = t**beta*np.exp(-K2*t**delta)
    g2 = -K1*t**(delta + beta)*np.exp(-K2*t**delta)
    g3 = K1*t**beta*np.exp(-K2*t**delta)*np.log(t)
    g4 = -K1*K2*t**(delta+beta)*np.exp(-K2*t**delta)*np.log(t)
    return np.array([g1, g2, g3, g4]).T


fig, ax1 = plt.subplots()
popt, pcov = curve_fit(FFA, t, y, jac=grad, maxfev=2000)


K1, K2, beta, delta = popt
print("value of K1 = ",K1)
print("value of K2 = ",K2)
print("value of beta = ",beta)
print("value of delta = ",delta)


ax1.plot(t, y, "o")
teval = np.linspace(t.min(), t.max(), 201)
yeval = FFA(teval,*popt)
ax1.plot(teval, yeval, color='red')

ax1.set_xlabel('Time (in s)')
ax1.set_ylabel('Spectral flux density (in Jansky)')
ax1.set(title='Light Curve')
plt.show()

结果确实表明您有一些缩放问题

value of K1 =  858036532.7666308
value of K2 =  21.214387300160098
value of beta =  5.958499311376985
value of delta =  0.3661586949061432

我建议您事先对模型进行无量纲化，以确保您的所有数字都处于相同的数量级。

更新：2021 年 11 月 12 日

您更改了模型，但我将其重写为

f = K_{1} t^{β} e^{- K_{2} t^{- 3}} .

$f = K_1 t^\beta e^{-K_2 t^{-3}}\, .$

$(1.39/5)^\alpha$ 和 $(1.39/5)^{-2.1}$ 是固定的数字，可以被吸收到 $K_1$ 和 $K_2$ . 使用这个模型对我有用。

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

t = np.array([33.90, 76.95, 166.65, 302.15, 330.11, 429.82, 533.59, 638.19, 747.94])
y = np.array([0.25, 1.81, 8.32, 11.60, 12.18, 10.12, 9.44, 5.81, 5.42])
yerr = np.array([0.09, 0.08, 0.14, 0.13, 0.16, 0.11, 0.06, 0.05, 0.06])


def FFA(t, K1, K2, beta):
        return K1 * t**beta * np.exp(-K2 * t**-3)


fig, ax1 = plt.subplots()
popt, pcov = curve_fit(FFA, t, y, sigma=yerr, maxfev=2000)


K1, K2, beta = popt
print("value of K1 = ",K1)
print("value of K2 = ",K2)
print("value of beta = ",beta)


ax1.plot(t, y, "o")
teval = np.linspace(t.min(), t.max(), 201)
yeval = FFA(teval,*popt)
ax1.plot(teval, yeval, color='red')

ax1.set_xlabel('Time (in s)')
ax1.set_ylabel('Spectral flux density (in Jansky)')
ax1.set(title='Light Curve')
plt.show()

以下是结果

value of K1 =  13773.167296442545
value of K2 =  6542145.117006296
value of beta =  -1.1791460005485805

协方差矩阵是

[[ 4.90022065e+08  3.82085117e+10 -5.62753836e+03]
 [ 3.82085117e+10  4.62424755e+12 -4.33603182e+05]
 [-5.62753836e+03 -4.33603182e+05  6.47333340e-02]]

其它你可能感兴趣的问题

上一篇Ubuntu 上的 COMSOL（10.04 LTS 和 11.04）下一篇矩阵是如何存储在 MATLAB 中的？