如何將 SIR 和 SEIR 模型擬合到流行病學數據？
我有2014年塞拉利昂埃博拉病毒爆發的流行病學數據[ xlsx ] [ csv ]。
我想用SIR隔間模型來模擬爆發
和SEIR隔室模型
很像（聲稱是）在本文印刷版的第 7 頁上完成的。
由於 S(E)IR 模型不承認緊密的分析解決方案，我在 Matlab 中對它們進行了建模，每個模型都有幾個*.m*文件。
SIRandSEIR函數返回系統的累積狀態（我將其解釋為每個狀態變量的積分），因為我只有流行病學公告中的累積數據 。
先生
SIR_factory.m
% Returns a function that returns the evaluation of the ODE defining the % system at time t. % % Using a High-order function saves me from learning how to pass % parameters to the function used by the ODE solvers :P function ode = SIR_factory(beta, gamma, N) %Returns an array with [dS/dt dI/dt dR/dt] ode = @(t,y) [ %y(1) = S, y(2) = I, y(3) = R -beta*y(1)*y(2)/N beta*y(1)*y(2)/N-gamma*y(2) gamma*y(2) ];先生
% For some reason the function must handle the input x in a matrix form or % the fitting app invocation of this function will throw a matrix concat % error % %INPUT: % beta and gamma are the system parameters % S0 and I0 are the initial conditions % x is the time at which the state of the system is requested % i is the index of the state of the system to return (1 = S, 2 = I, 3 = % R) function R = SIR(beta, gamma, S0, I0, x, i) %Total population is S0 (Susceptible individual) plus I0 (Infected %ones), this means R0 is implicitly 0 N = S0 + I0; %Find the max of the matrix M = max(max(x)); %Pre-allocate the output R = zeros(size(x, 1), size(x, 2)); %Solve the system, from time 0 to the max time requested, with the %initial conditions given [t,y] = ode45(SIR_factory(beta, gamma, N), [0 M], [S0, I0, 0]); %Comput the "integral" of the state variable requested y = cumtrapz(t, y(:,i)); %Fill the output (Ugly) for r=1:size(R,1) for c=1:size(R,2) R(r,c) = interp1q(t, y, x(r, c)); end end
詩
SEIR_factory.m
% Returns a function that returns the evaluation of the ODE defining the % system at time t. % % Using a High-order function saves me from learning how to pass % parameters to the function used by the ODE solvers :P function ode = SEIR_factory(beta, gamma, delta, N) %Returns an array with [dS/dt dE/dt dI/dt dR/dt] ode = @(t,y) [ %y(1) = S, y(2) = E, y(3) = I, y(4) = R -beta*y(1)*y(3)/N beta*y(1)*y(3)/N-delta*y(2) delta*y(2)-gamma*y(3) gamma*y(3) ];SEIR.m
% For some reason the function must handle the input x in a matrix form or % the fitting app invocation of this function will throw a matrix concat % error % %INPUT: % beta, gamma and delta are the system parameters % S0 and E0 are the initial conditions % x is the time at which the the state of the system is requested % i is the index of the state of the system to return (1 = S, 2 = E, 3 = % I, 4 = R) function R = SEIR(beta, gamma, delta, S0, E0, x, i) % The total population is the number of Susceptible individualts plus % the Exposed one. This implicitly sets I0 = R0 = 0 N = S0 + E0; %Max of the input matrix M = max(max(x)); %Pre-allocate the output R = zeros(size(x, 1), size(x, 2)); %Solve the system from time 0 to time M with the initial given %conditions [t,y] = ode45(SEIR_factory(beta, gamma, delta, N), [0 M], [S0, E0, 0, 0]); %Compute the "integral" of the required state y = cumtrapz(t, y(:,i)); %Fill the output matrix for r=1:size(R,1) for c=1:size(R,2) R(r,c) = interp1q(t,y,x(r, c)); end end目的是使這些模型適合爆發的數據。
為此，我將報告日期轉換為自第一份報告以來的天數，並為每一行計算感染人數，即Total cases減去Total Deaths。
生成的 CSV 在這裡。
作為*.xlsx*文件的原始數據在此處。
為了了解數據，這裡是摘錄
... 121,235 125,301 128,300 131,322 132,373 136,419 140,415 141,449 143,462 147,483 149,533 150,518 156,604 165,770 ...在添加適當的路徑後，我將這些數據讀入 Matlab，如
>>>DATA = csvread('ebola-sierra-leone-inf.csv') >>>DAYS = DATA(:,1) >>>INF = DATA(:,2)
然後我使用曲線擬合應用程序來擬合模型。
擬合會話*.sfit*文件在這裡。
在 2014 年，塞拉利昂人口約為 707.9 萬，我將其用作初始減去代表的未知量(先生) 或(SEIR) 變量。
先生
我創建了一個SEIR擬合，使用X
DAYS數據和Y數據。INF我選擇了一個自定義表達式方程，
SIR(b, c, 7079000 - I, I, x, 2)它代表一個帶有參數的 SEIR 模型b() 和c() （**全部限制在 [0, 1]**中，初始值為 0.5）並返回隨著時間的推移 感染個體的累積數量。擬合參數
I被限制在 [0,)，初始值為 1。詩
我創建了一個SEIR擬合，使用X
DAYS數據和Y數據。INF我選擇了一個自定義表達式方程，
SEIR(b, c, d, 7079000 - E, E, x, 3)它代表一個帶有參數的 SEIR 模型b(),c() 和d() （**全部限制在 [0, 1]**中，初始值為 0.5）並返回隨時間推移 的累計感染個體數。擬合參數
E被限制在 [0,)，初始值為 1。
結果
Matlab 發現非常不合適。
先生
General model: f(x) = SIR(b, c, 7079000 - I, I, x, 2) Coefficients (with 95% confidence bounds): I = 21 (16.63, 25.38) b = 0.9999 (-0.1879, 2.188) c = 1 (-0.1859, 2.186) Goodness of fit: SSE: 3.124e+08 R-square: 0.9153 Adjusted R-square: 0.9147 RMSE: 1092
詩
General model: f(x) = SEIR(b, c, d, 7079000 - E, E, x, 3) Coefficients (with 95% confidence bounds): E = 43.96 (-20.21, 108.1) b = 0.9989 (-2.32, 4.318) c = 1 (-2.317, 4.317) d = 1 (fixed at bound) Goodness of fit: SSE: 3.252e+08 R-square: 0.9119 Adjusted R-square: 0.9112 RMSE: 1114
然而…
…如果我將人口數量從大約 707.9 萬更改為 70.790（原始的 1/100），兩者的擬合都會提高
顯然，我在概念上和/或數字上做錯了什麼。
如何將 S(E)IR 模型擬合到流行病學數據？
我將把我的評論限制在 SEIR 模型上——SIR 模型的問題是相似的，無論如何它都可以被視為 SEIR 模型的一個特殊限制情況（對於大型).
到目前為止你做了什麼
我看過你的 MATLAB 代碼，這對我來說絕對沒問題。對於給定的一組模型參數，您的代碼求解 SEIR 微分方程以給出函數,,,在某個時間間隔內。然後計算累積狀態它用作擬合模型的基礎（如果我在這裡錯了，請糾正我）。
可用數據：你有時間序列和，分別是累計病例數和死亡數。通過最小化曲線之間的差異來進行模型擬合和. （這假設一個疾病病例對應於一個過渡到狀態。）獲得了較差的擬合。置信區間的意義也值得懷疑——即使所有模型參數都被限制為正，下限通常為負。
模型與數據
我可以看到指定的 SEIR 模型與可用數據相關的方式存在幾個問題。第一以上不代表感染人數，簡單來說. 看來您實際上想要等同於到. 計算似乎沒有必要。
其次，您似乎隱含地假設“恢復”（過渡到類別）總是導致死亡。然而，我理解——就埃博拉而言——也有可能被“治愈”。因此，可用的死亡數據不能與您設置的 SEIR 模型中的變量直接相關。這表明需要一種模型來考慮埃博拉病毒可能出現的不同恢復模式。
第三個問題是，通過從另一個數據時間序列中減去一個數據時間序列，您會丟失原始數據中的一些信息。理想情況下，最好使用兩個可用的時間序列來擬合模型。
修改後的 SEIR 模型和擬合程序
為了改進模型擬合，我建議查看本文中所做的建模。在這裡，他們為埃博拉病毒使用了修改後的 SEIR 模型，看起來像
這裡是病死率，所以狀態對應於“治愈”。在該模型的上下文中，累計病例數為累計死亡人數為. 也許可以在 MATLAB 中同時擬合這兩條曲線？ 其他型號
更複雜的模型當然是可能的，例如參見本文，其中使用了其他疾病類別。我們還可以添加隨機性、更詳細的接觸結構模型等。將傳播模型擬合到 2014 年埃博拉疫情數據是一個活躍的研究領域。不過，您可能希望使用上面修改後的 SEIR 模型獲得合理的擬合。我想說的是，將傳播模型擬合到埃博拉疫情數據不是一項簡單的任務！
最後：您提到的論文似乎不是同行評審的期刊文章。它也是匿名的。我不會依賴它作為信息來源。