" O' z9 L4 V0 U$ f# C& z& h, e数学建模之传染病SIR模型(新冠真实数据)) z7 M% {# w# ]+ T
传染病模型的基本问题 8 A# g* ^; g9 |: L6 l4 ^. E" t ~4 I描述传染病的传播过程# w- p: q4 d$ C, c5 z
分析受感染人数的变化规律 ! e! i, m; Q3 y6 f1 h& j1 Q* A预报传染病高潮到来的时刻! t! H! h! p* O5 ~& z8 @
预防传染病蔓延的手段. y) s, [! g/ Z" _# M8 b4 d$ z) O
按照传播过程的一般规律用机理分析方法建立模型2 {2 K0 P$ J5 j
注:我们这里是介绍数学医学领域中基本的传染病模型。不从医学角度分析各种传染病的特殊机理,按照传播过程的规律建立微分方程模型.# l4 K, l- c6 Y9 I
. p% Y3 H) @: k8 H. h' T1 S 2 T8 m" M: K4 q- M- `建立模型; O1 Y# t* W" B8 A( ~6 J
模型一 , C; _" H% i# V; ^% v假设: 7 V( A T. X8 M & \1 u8 b1 Y6 x+ R# c" D* o+ \/ _0 L3 c2 V, p: R S, [6 c0 w
设已知感染人数为i ( t ) i(t)i(t)(病人数量随时间变化) - t& c8 G+ P& g5 o, [, i设每个病人(单位时间)每天有效接触(足以使人治病)人数为λ \lambdaλ 8 q5 f) A# N6 ~" M$ D- k) I模型:; q% r, ~' X8 r \# {$ r3 E
单位时间Δ t \Delta{t}Δt内,新 增 的 人 数 ( 现 有 − 原 有 ) = 原 有 的 × λ 新增的人数(现有-原有)=原有的 \times \lambda新增的人数(现有−原有)=原有的×λ,即1 n1 b7 I( _- d. z1 e( L: x
$ x+ M2 G& K: T : V Z7 ~2 k Z' ri ( t + Δ t ) − i ( t ) = λ i ( t ) Δ t i(t+\Delta{t})-i(t)=\lambda i(t)\Delta{t}i(t+Δt)−i(t)=λi(t)Δt7 A3 I' s3 t/ L+ @: S! B
一开始的感染人数为i 0 i_0i 7 y; \" x z$ Y7 y( Q% ?( e8 r2 f
02 w6 Q! r0 U M( o
% e+ Y; p! R4 w, {
; a% b6 I2 U( B& u1 n
i ( 0 ) = i 0 i(0)=i_0i(0)=i 0 G$ o3 p9 Z0 v$ F( V7 F
0. P8 t3 F1 x. K' Y) m+ q
T* E- z( n5 ~ z( a4 T7 l, m3 b
/ t3 b+ X8 [$ U: z4 |0 e8 P" S解微分方程可以得到 7 K; [. Z+ r9 H) ~+ I. F- z7 G& gi ( t ) = i 0 e λ t i(t)=i_0e^{\lambda t}i(t)=i + `) l' m a! h% f- b! c/ G
0 c! z9 l2 Z9 G/ g" n1 Z* H 2 _8 A4 a+ b6 z; a e ; x* \6 N8 g) `- a) ^" s @7 K7 H) Hλt/ g& I3 e7 m/ w% e
0 s/ B- t9 y9 R$ O9 \! k; W8 X+ Z( W
所以可以可到当λ → ∞ \lambda \rightarrow \infinλ→∞时i ( t ) → ∞ i(t) \rightarrow \infini(t)→∞ 1 J* _% _# r+ I2 S% S+ E当然这是不可能的,因为我们考虑的因素太少了,首先一个是,若有效接触的是病人,则不能使病人数增加,所以必须区分已感染者(病人)和未感染者(健康人)看模型二来解决这个问题8 I4 v$ F& F4 t; Z- ]
4 S. U, u* m. O2 E2 ^' [
+ \- f4 e5 [9 r, ~ \8 `
模型二4 }5 ~0 ^* j9 y1 Z5 A
假设:( x4 c' @2 c0 _- g8 r, ]
1 r7 F5 F# @* R5 h1 o8 _
$ \0 @2 L) z3 M6 o& v. Q4 k将人群分为两类:易感染者(Susceptible,健康人)和已感染者(Infective, 病人).7 v( e- V+ B3 @
总人数N不变,时刻t健康人和病人所占比例分别为s ( t ) s(t)s(t)和i ( t ) i(t)i(t), 有s ( t ) + i ( t ) = 1 s(t)+i(t)=1s(t)+i(t)=1 0 W0 p! I6 U3 O( z! L每个病人每天有效接触人数为λ \lambdaλ(日接触率),且使接触的健康人致病.* G6 A h4 |% q4 B
建模:% X7 Z$ Y6 H4 G- K* e
每天新增的总人数为原有的人数乘以每个人可以传染的健康的人数,再乘Δ t \Delta tΔt o% Q& f, r* _, l ? c7 n0 h/ `" @: ^9 b* L
1 q3 }+ q; d/ v: f" b. H; e
Δ t \Delta tΔt除过去,两遍N约分得到下面,( R3 J+ B+ P' g: W% |! }
+ d6 Z. k* }& y/ N2 W1 y6 a4 E# m ' r3 A) m. e. U, f( Z& h: \( SMATLAB解一下这个微分方程0 }3 T# f, k! u! }
* x. ?8 T7 n ^# e, U9 ]
4 `) h( n) }" j( ey=dsolve('Dy=n*y*(1-y)','t'); : H, M2 R* @, i- D1 E , f2 i' @+ L% W. C' \1 l, h2 } ( {+ k, r, `: Ay = 0 H N4 F. i0 M% @$ [5 G/ R: X# M -1/(exp(C1 - n*t) - 1)/ s9 E9 f7 k5 h* P$ T+ U6 q
0. \( S& O5 m @- |2 y+ r3 b$ j
1 6 r4 O7 S! ^% Z1. c8 x6 c+ \6 O4 b, z2 U0 o/ J
2) C) \0 U) x+ Y [/ k
3 ' x7 I' ^" k- k; d7 J4 . N; ?$ y* T1 n; O n5, o0 ^4 H y. g9 g, r
6 l" o+ Y1 S2 s. j写规范点就是这个函数# m8 c* v8 F7 }) G
0 I" J* Z: L e9 q: U' Z
" a3 T* L& a- M1 \/ U函数图像大致为8 J/ S4 V! p, h! T0 O
+ Q0 e s+ y% s% @" o" O1 R3 }; x) f2 B % J8 u% j/ V; O: R' B% h) v可以看出t = t m t=t_mt=t * }$ d. R8 J0 e( zm ( N# _& e, z: l" N$ P' Z - [2 ^4 M- ? X: z5 N 时这里图像的斜率有个最大值,其也就是传染的最快的时候,即传染病的高潮时刻,当然t m t_mt " k' j3 s% t9 `. `
m . l( Q" \( b) v) G ' [* E L$ J# J6 ? 是可以求出来的 + v. L* |, Z! _& _6 p9 Y0 `' z' Q1 U
( x9 Y) b. o0 F' K再看原式,当t → ∞ t\rightarrow \infint→∞时i → 1 i\rightarrow 1i→1# V) _1 E1 L# ?! D7 ?, C
病人的比例为1,当然这也是不可能的,因为我们还没有考虑有没有可能治愈,看模型三 8 l# _, D. ^/ z: r( B! @9 n" j' Z6 s" y. |. P2 S
+ b6 ]1 o2 ~: \ & ~( h T5 I) p实战建模 ' m$ e6 S' \3 r4 X$ Y数据处理% I. b7 r3 x% V7 `
5 Z7 M7 w9 M5 v9 S* g& A
' y- W- x4 I* Y9 O4 `. m首先,我用python爬虫爬取了丁香医生官方数据,一共5534条数据 特征包括感染、死亡、治愈的总数,当日感染、死亡、治愈新增,疑似病例,时间,省份等14个特征4 d/ I) B/ r( }9 D1 m" l+ l
: t, V2 s) \0 E0 D# ?& [0 K
. c. g2 j) D: V. O/ q/ [# {
% B" s9 G+ e* q) U% a, L7 N, q" U! ?6 J" h3 d5 L
然后用python进行数据提取,提取了较为典型的湖北省的数据作为我的参考依据 . ~& I) z7 h5 r3 C" m5 ?, {+ B" ]/ M6 {- T$ k6 |
0 V$ T4 j" k i# R4 S
$ G5 m/ X. A, t0 e& ]: s( H4 J8 y: a" z$ p; h
然后用python对数据进行清洗,提取出了患病总数,现存患者总数,死亡总数,治愈总数,时间,省份这几个特征3 ~) w$ S" @' j1 g5 S5 N+ c
0 b; s' R1 q* z5 N+ \9 y8 t3 p4 I6 a6 a
对日期格式进行修改,值保留月和日,并与死亡人数的位置交换( r0 g) _8 K% f1 E
6 g* x+ N- E c. X 9 f4 f# b/ @2 n9 u1 v ]# e这里我用python对提取的四个特征分别进行了数据分析(主要包括计算最值,平均值等,),并把1.20日作为第一天,7.02日作为最后一天也就是第165天,做了可视化可视化处理。. |* o9 U7 f5 B* H
感染人数示意图; D: u6 ?3 L3 f. ?, e
/ |5 E% U9 Z& E2 D S
- Y3 T+ }& Y) |. x# J% f
治愈人数示意图: T" u; @; C" m, a- \' ~
6 \. @& ]; W1 u) C8 E8 F
" y9 n7 ^, \% Z: K! }! V- [
% b% h4 @# H- J: ~2 L3 v; ]7 y3 R( q) f4 ]# ~& G
现存患者数量图 9 D! ?1 r M. @6 B' X , `8 O& R1 T& J- \( V; Y3 U6 G* t; ^/ m) M. q) V4 A( j
死亡人数示意图2 h4 C* w c4 J; |1 [" C7 K
. h+ K# L" S, n& u7 @5 [# l
" H. P% n) N. `
# J$ [% D& I% C# o
" ^" c$ y6 I8 {经过上面的图片与describe数据分析,我们发现有一天是异常的,患者多出了平时的十倍左右,经过查阅资料,这天因加强了检测标准,所以增多了很多。为了避免这个数据的影响我们选择将这一天删去(或者用平均数或中位数代替也可) : |6 y+ z C: J P7 O将上面清理过的数据存放到csv文件中6 ?7 k" A( a! x3 h
. u5 e& F/ ~+ @+ ^4 Y8 q8 z! i7 N' R+ h r) g
模型建立( ^6 v7 C8 h' D& R
模型假设 % \/ A5 e! J/ y' y" [) z, a/ g经过上面数据的分析,我们大体可以进行如下假设:7 V' e$ ]1 s* w0 C$ \- ?
1.由于不存在封闭情况,考虑开放体系。 . | Q+ k' y1 x3 E' R2.目前数据以天为单位发布,因此不考虑连续变化情况,只考虑离散的方程。; @$ \) A: b5 H% [$ j+ {6 `- w
3.新型冠状病毒的治愈人数和死亡人数相对较 小,因此只考虑 Susceptible(易感)和 Infected(感染) 两类人群。设易感人群总数为N 9 j9 _8 F% x' d* |0 [4.经专家鉴定新冠病毒患者治愈后至少六个月之内不会再被感染,所以设治愈后移出易感人群。4 W2 t$ f1 M9 A7 `& i1 k
5.设每个病人每天有效接触人数为 λ \lambdaλ(日接触率),且使接触的健康人致病. ' t, z: r% k0 b) H( ?6.设病人每天治愈的比例为 μ \muμ(日治愈率). \3 g. X1 F8 r( v* n1 o t4 m" @% O
7.时刻t健康人、病人和移出者的数量分别为 s(t), i(t), r(t). W! {% H1 M* n; l5 [
" y, l& ]* G! F4 b5 K; R: {* H9 N0 X6 [! N1 j( n9 h7 Z* Q
模型一. X4 K0 O: F% _5 x: i" ]# F3 e
+ k1 ^, k+ u D( {/ n
g0 ?+ @( B7 F1 S |( i分析可以得到移出者r(t)=治愈人数+死亡人数5 s! a4 x( W; N3 V/ X) A
通过python数据处理,我们算出了r(t)的值,并将其可视化 ; g: }! g8 O5 Z6 L, C o, {7 L( [, ^+ x3 g& V3 H! M; z! Y) a' W
! A1 |# h5 _- A( e2 n2 o实际上,λ \lambdaλ (日接触率)和 μ \muμ(日治愈率)都是随着时间变化的,这里我们设s(t), i(t), r(t) 为第t天健康人、病人、移除者(病愈与死亡之和)的数量, s(t)+ i(t)+r(t)=N.. 5 ]9 V5 c9 U! n0 X5 B r p(t), (t) ~第t天感染率, 移除率(治愈率与死亡率之和) ) {* M: O( D. I" ?( Y) X* O/ j2 R0 R有 d i / d t = λ ( t ) s ( t ) i ( t ) − ( t ) i ( t ) di/dt=\lambda (t)s(t)i(t) - (t)i(t)di/dt=λ(t)s(t)i(t)−(t)i(t) 3 w% @6 n+ e, S' U4 I( D7 M" v因为s远大于i, r,s(t)视为常数,所以有) m% |+ q/ E n/ P+ N" i6 l