时间戳

首页 > 解决方案 > 优化/曲线拟合

python - 优化/曲线拟合

问题描述

我正在为 covid 数据拟合流行病学模型。代码如下所示。该模型有 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,95.9,111.4,143.7,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, beta, 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, t, args=(N, beta, gamma))
S, I, R = ret.T
plt.plot(t, I, color='red', linewidth=2)
# plt.savefig('fitted.pdf', dpi=300, bbox_inches='tight')
error = sum(I - c)

标签： pythonoptimizationcurve-fitting

解决方案

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

推荐阅读