优化/曲线拟合

如何解决优化/曲线拟合

我正在采用流行病学模型来收集数据。代码如下所示。该模型具有3个参数：N，beta和gamma。我如何优化这些参数以得到最佳模型拟合（按“错误”衡量）？

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

# Data
t = np.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158])
c = np.array([111.0,138.3,135.3,143.7,132.0,114.1,95.9,76.4,47.9,56.9,53.6,51.4,56.6,66.4,67.4,71.1,63.3,74.0,83.7,95.1,94.4,111.4,140.7,146.9,150.0,161.3,192.6,205.3,189.6,198.3,213.3,218.7,245.0,242.0,247.3,282.1,319.6,346.7,368.0,351.1,369.3,420.6,440.6,461.7,490.3,539.7,609.4,643.9,661.3,705.0,785.7,825.9,835.7,847.0,914.0,943.0,998.3,1009.7,1026.0,1009.1,1133.4,1180.9,1302.1,1374.9,1442.9,1534.6,1649.7,1655.4,1912.7,2027.7,2143.7,2228.9,2360.3,2454.1,2556.6,2539.4,2641.9,2823.3,3107.6,3236.3,3334.7,3570.1,3617.4,3684.0,3849.3,3894.9,4008.1,4253.4,4483.4,4926.4,5267.9,5588.4,5833.1,6096.3,6443.0,6791.0,7098.0,7504.9,8025.3,8373.7,8779.6,9235.1,9333.1,10039.7,10509.0,10886.7,11356.0,11725.0,11776.7,12340.6,12268.9,12415.3,12385.0,12583.7,12261.7,11929.4,11985.6,11975.9,12057.4,11903.0,11586.4,11271.6,11137.6,10882.1,10588.1,10169.6,9870.0,9436.0,9190.4,8793.9,8393.4,8002.1,7470.4,7128.3,6910.6,6676.6,6398.7,5577.4,4954.4,4809.1,4352.1,3926.6,3755.4,3719.3,3877.3,3867.9,3456.9,3341.7,3204.0,3080.6,2981.9,2805.9,2620.9,2399.1,2215.1,2183.3])
plt.title(r'Daily cases - SIR fit: N=42 000,$\beta=0.148$,$\gamma=0.05$')
plt.xlabel('Days')
plt.ylabel('Cases (7-d moving average)')
plt.plot(t,c,".")

# Model
# Total population,N.
N = 42000
# Initial number of infected and recovered individuals,I0 and R0.
I0,R0 = 1,0
# Everyone else,S0,is susceptible to infection initially.
S0 = N - I0 - R0
# Contact rate,beta,and mean recovery rate,gamma (in 1/days).
beta,gamma = .148,1/20
# The SIR model differential equations.
def deriv(y,t,N,gamma):
    S,I,R = y
    dSdt = -beta * S * I / N
    dIdt = beta * S * I / N - gamma * I
    dRdt = gamma * I
    return dSdt,dIdt,dRdt
# Initial conditions vector
y0 = S0,I0,R0
# Integrate the SIR equations over time t.
ret = odeint(deriv,y0,args=(N,gamma))
S,R = ret.T
plt.plot(t,color='red',linewidth=2)
# plt.savefig('fitted.pdf',dpi=300,bbox_inches='tight')
error = sum(I - c)

解决方法

from scipy.optimize import minimize

def fn(x):
    # parameters unwrapped
    N,beta,gamma = x;
    
    # Initial number of infected and recovered individuals,I0 and R0.
    I0,R0 = 1,0
    # Everyone else,S0,is susceptible to infection initially.
    S0 = N - I0 - R0
    # Initial conditions vector
    y0 = S0,I0,R0
    # Integrate the SIR equations over time t.
    ret = odeint(deriv,y0,t,args=(N,gamma))
    S,I,R = ret.T

    # absolute error to avoid negatives
    error = sum(abs(I - c))
    
    return error

# initialise with current best guess
init_x = [42000,.148,1/20]

# calculate result
res = minimize(fn,init_x,method='Nelder-Mead',tol=1e-6)

# calculate final error
fn(res.x)

由于您具有自定义函数，因此最小化此函数的一种简单方法是调用scipy's minimize参数。该函数仅接受一个参数，因此这三个参数需要作为列表传递。请参阅上面的示例。然后，您可以通过调用res.x来获得最佳参数。

我建议您看看Scipy's optimize.curve_fit函数。在这里，您可以定义可以适合您的数据的功能。返回最佳参数和拟合协方差。

如何解决优化/曲线拟合

解决方法

相关推荐