新冠疫情預測

2020年2月春節先後, 新冠狀病毒疫情打亂了全部中國人對春節假期的安排和對新年的嚮往, 我也是。原本是計劃去東北的大城市走一趟的。 一家人都悶在家裏一個月，沒人知道疫情什麼時候結束，何時能重新開始工做，開始上學和正常的生活，全部人都 很焦慮不安。因而向用一個有效的傳染病疫情模型預測一下。

 1.　　SIR模型簡介

新冠狀病毒疫情同2003年的非典疫情相似， 基本上從人羣劃分來看， 有以下特色：

• 健康人沒有免疫力所以可能會被感染病毒，可稱之爲易感者
• 感染者具備傳染性，所以能夠傳染給健康人將其變成感染者。
• 感染者或者死亡從而將被活化掩埋處理，或者治癒從而得到免疫力。總之死亡者和治癒者再也不可能會感染病毒， 可稱之爲疫情移除者 。

可用SIR （Suspectible，Infective and Removed )這三個字母來表明以上的三我的羣： 易感者，感染者，疫情移除者。SIR 模型的目的就是對三類人羣數量的變化創建數學模型，以此來預測疫情什麼時候到達拐點 （感染者每日淨增量由遞增變爲遞減），峯值（累計感染者達到最大值或天天的感染者增量爲零）和結束（累計感染者所有轉化爲移除者，淨值爲0）。對照衛健委天天發佈的官方數據，以下所示， 包括確診（Infected), 疑似（Suspected), Dead(死亡）, Recovered(治癒），Observed (醫學觀察），Contacted (密切接觸）, SIR 的源數據以下：

• Removed （移除者）是累計死亡和治癒之和：Dead + Recovered
• Infective（感染者) 是累計確診數減去移除者數 : Infected - Removed 
• Suceptible (易感者）： 全部的除感染者和移除者以外的中國人都屬於易感人羣 ：人口總量 N（14億）- Removed - Infective 
• Susceptible(S) + Infective (I) + Removed (R) = N

SIR 模型假設：

1. 人口總數（N）保持不變，N 由 S易感者人數(S），感染者人數（I) 和移除者人數（R） 組成，N = S + I + R 
2. 易感者(S)與感染者(l）有效接觸後，天天按必定的比率從易感者成爲新的感染者。這個比率稱爲日接觸率λ （讀lambda) .
3. 感染者通過治療後天天按照必定的比率μ(讀miu）從感染者成爲移除者（R康復或死亡）， 這個比率稱做日治癒率或日移除率。
4. 移除者（Ｒ）不會再次感染病毒或傳染病毒給別人。
5. 感染者（I）天天的增量是易感者(S)天天轉化成I 與 天天I 轉化成R的淨值。

注意以上描述都是在圍觀描述Ｓ，Ｉ，Ｒ天天的增量或減量，在數學上能夠用你微分方程來表示，以下所示：

dS/dt = -λ ·S/N · I          (1)

dI/dt= λ ·S/N · I - μ·I       (2)

dR/dt = μ ·I                     (3)

N = S + I + R                  (4)

雖然數學語言高度抽象，但以上方程和5點假設描述的是一樣的意思。好比衛健委的數據是以天爲單位發佈的， 因此dt = 1(天）是天天的意思。

dR/ｄt 表示移除者R天天的增減：好比移除者2月1日相對於1月31日的增量，這個增量與當天的I成正比，比率爲日治癒率μ. 

dS/dt， 表示易感者天天的增量，它是一個負數，數值與當前易感者在總人口占比和當前感染者比正比，比率爲日接觸率λ 。

dI/dr , 表示感染者天天的淨增量， 他是dS/dt的絕對值 和 dR/dt 差

S， I ， R　各自的數值天天都不同，但三者之和爲人口總數N。N是常數。

此處最好有個圖來描述，讓這個變化過程形象起來，讀者本身想象吧。

 S，I，R都是關於ｔ的未知的連續函數，若是能求出他們各自的解析解（S，I, R 的各自函數形式），而後將他們的圖形畫出來，那麼將來趨勢就一目瞭然了吧啊，不幸的是解析解沒法求出，所幸的是咱們的目的是預測將來一年的趨勢，數值解一樣能達到目的。

 數值解是經過給定初值T_０時Ｓ_０，Ｉ_０，Ｒ_０，並利用微分方程組的微分等式，遞歸的計算將來（好比２０２０年１月到１２月）每一天的Ｓt，Ｉｔ，Ｒｔ　．

那麼接下來的問題是肯定微分方程組的參數 λ和μ了。

2.　常數接觸率和治癒率的模型實現

ＳＩＲ模型假設接觸率λ和治癒率μ是常數的 , 天天Ｓ, I, R 的數值和他們的微分的可計算以下 .

R =  累計治癒+ 累計死亡, I = 累計確診 - R, S = N - I - R.

ｄR/dt = R_t+1 - R_t, ｄI/dt = I_t+1 - I_t , ｄS/dt = S_t+1 - S_{t .}

 

 

 而後將S, I, R, dR , dI 帶入公式(2)和(3) , 利用最小二乘的優化方法可求出最優的常數接觸率λ和治癒率μ . 而後在求出將來一年的常微分方程組(ODE）的數值解，再將這些數值在日期軸上畫出來，就能開到將來了．下面是部分的matlab代碼, 具體詳見github .

 

 sir_run.m

load('weijianwei');   %load the source data 
u = regress(dR(1:end-1), I(1:end-1) ); % optimize u
lamda = regress(dI(1:end-1)+ dR(1:end-1), S(1:end-1).*I(1:end-1)/N); % optimize lamda

%lamda=0.492112631178384; %optimized by using 1stopt
%u=0.421005009555275; %optimized by using 1stopt

%numeric solution of SIR ODE 
[t,y]=ode45('sir', [0,400], [1399999802    170    28]);

% plot the the full SIR graph 
...

sir.m : sir ODE 

%S' = -lamda * S * I / N , susceptible
%I' = lamda * S * I /N - u * I , infetive , in hospitable 
%R  = u * I  , recovered or dead , not in hospital 
%S  +I + R = N

function y=sir(t,x)
    global lamda u N
    y=[-lamda*x(1)*x(2)/N,lamda*x(1)*x(2)/N-u*x(2), u*x(2)]';
    

　　運行一下看看效果

 

 優化的效果很不理想, 14億人口在5月前所有感染,高峯期感染11億, ----- 這個不可能的 .....

看一下具體的日接觸率λ和日治癒率μ 是多少呢, 以下

 

日接觸率遠大於治癒率,大概是24倍 , 這樣的比率會形成易感者迅速大量感染，因爲易感人羣趨於零，感染者人數達到頂峯,　而後慢慢治癒 , 者就是上圖所顯示的樣子  .

在看看數據擬合的程度怎麼呢?

 

 

 

 

 

 

 雖然沒有作擬合率的量化計算(好比R，假設檢驗等），目測擬合度是很是低的．　

看起來假設lamda 和u爲常數是行不通的 ，想別的辦法吧 。

 

3.　帶衰減的接觸率和治癒率的模型實現

 lamda 和u爲常數的假設通過截止到2月2日新冠狀病毒疫情的數據驗證，是失敗的，聽說在２００３年的 非典疫情，這個假設也是失敗的。那麼lamda 和u也是關於ｔ的函數嗎？　根據SIR 常微分方程的（２）（３），　能夠用一下公式求出天天的日接觸率λ和日治癒率μ，看看他們是否是顯得隨着時間有規律的變化．

u_t = dR ./ dI 
lamda_t = (dI + dR) ./ (I .* S / N)

下圖黃色線是lamda_t,  藍色線是有指數函數:a·e^b·t , a 和  b 是常數，t 表明從1月19日到改天的天數。雖然擬合也不是很好,感受上27以前的噪音較大, 27以後有規則的指數形式衰減。儘管受噪音影響，圖中的lamda decay 曲線仍是比較不錯的擬合了lamda隨時時間而逐漸衰減的現實。 不過感受衰減的過快，十幾天lamda 衰減爲初值的一半，衰減過快會使Infective 曲線機制過早到達。

 

 一樣的， u_t治癒率 也有隨時間以某種指數形式遞增的趨勢， 但感受遞增的有些快。

 

 

逐漸衰減的lamda_t和逐漸上升的u_t, 會使逐步降日低接觸率和日治癒率的差距， 減小感染人羣的增加趨勢， 下降易感人羣的感染機率， 下降了I曲線的曲度，從而使極值下降。

下面是部分實現代碼

sir_t_run.m : 

clear;
global a b  N % a,b are the coefficients of decay function of lamda and u
N = 1.4e+9 ;
load('weijianwei');
u_t = dR ./ dI ; % u(t)
lamda_t = (dI + dR) ./ (I .* S / N); %lamda(t)

%plot lamda(t) and u(t)
dates = datetime(Date);
figure;
plot(dates, lamda_t, '--', dates, u_t, '--');
...

% fit lamda(t) using exponential 
x0 = [1, 0 ] ;
x = 1:1:length(lamda_t)-1;
a = lsqcurvefit(@decay,x0, x, lamda_t(1:end-1));

...

% fit u(t)
x0 = [1, 0 ] ;
x = 1:1:length(u_t)-1;
b = lsqcurvefit(@decay,x0, x, u_t(1:end-1));

...

% numeric solution of SIR ODE with decayed lamda and decayed u 
[t,x]=ode45('sir_t', [0,100], [1399999802	170	28]);
dates=start_date + t;

% plot the the full SIR graph 
figure;
plot( dates(1:60), x(1:60,2), '--', dates(1:60), x(1:60,3), '-.k') ;
...

sir_t.m : ode with decayed lamda and u 

%S' = -lamda * S * I / N , susceptible
%I' = lamda * S * I /N - u * I , infetive , in hospitable 
%R  = u * I  , recovered or dead , not in hospital 
%S  +I + R = N

function y=sir_t(t,x)
    global a b N
    lamda = a(1) * exp( t * a(2));
    u = b(1) * exp( t * b(2));
    y=[-lamda*x(1)*x(2)/N,lamda*x(1)*x(2)/N-u*x(2), u*x(2)]';

　　

運行一下看看效果吧

淨感染者在2月26左右到達峯值，極值爲6.5萬左右，加上移除者2.7萬左右的移除者，當時的累計確診應該在9.2萬人左右，疫情將在4月底完全結束。

 

 

 

 

 擬合的效果怎麼樣呢？還很不錯的，因此對上面的預測相對有些信心。 

 

 

 

 

 

 

4. 　SIR模型 變種

 衛健委的數據源裏還包括了一些重要的公開數據，如疑似病人，醫學觀察人數，密切接觸者，若是能見這個向量加入模型，他們的將來走向，也能反應疫情將來的發展。可是我不是很理解這些數據確切的含義，以及他們之間的關係，因此沒法正確的用ODE (常微分方程 ）對他們建模 。

好比疑似病人的去向是什麼，他來自醫院觀察者嗎？醫學觀察者是累計量仍是淨值？ 它的去向是什麼 ？ 它怎麼被選中的呢？

我在代碼裏試着創建一個模型（seir.m)除了SIR 還包含了疑似者（E), 醫學觀察者（O），結果不怎麼合理， 我相信模型是錯誤的。

但若是能知道他們的正確的定義和他們之間的轉化關係， 就能夠簡歷正確的微分方程， 從而創建更完善、更符合新冠狀病毒疫情的模型。

SIR仍是過於簡單，　好比SIR模型假設只有確診的感染者具備傳染性是不許確的， 由於處於潛伏期的疑似病人，醫學觀察者，密切接觸者均可能傳染性。若是能將他們加入模型，會更有助於疫情的預測。 

 

5.　　本文代碼

 https://github.com/kaixin1976/coronavirus