Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# coding: utf-8
2"""
3Implementation of a model for epidemics propagation.
4"""
5import numpy.random
6from sklearn.linear_model import ElasticNet
7from ._base_sir import BaseSIR
8from .covid_sird import CovidSIRD
11class CovidSIRDc(BaseSIR):
12 """
13 Inspiration `Modelling some COVID-19 data
14 <http://webpopix.org/covidix19.html>`_.
15 This model considers that observed data are not the
16 true ones.
18 .. math::
20 \\begin{array}{rcl}
21 S &=& S_{obs} (1 + a)\\\\
22 I &=& I_{obs} (1 + b)\\\\
23 R &=& R_{obs} (1 + c)\\\\
24 D &=& D_{obs}
25 \\end{array}
27 Where :math:`S`, :math:`I`, :math:`R` are
28 :math:`R_{obs}` are observed.
29 hidden, only :math:`S_{obs}`, :math:`I_{obs}`,
30 The following equality should be verified:
31 :math:`S + I + R + D = N = S_{obs} + I_{obs} + R_{obs} + D_{obs}`.
32 We also get from the previous equations:
34 .. math::
36 \\begin{array}{rcl}
37 dS &=& dS_{obs} (1 + a) = - \\beta \\frac{IS}{N} =
38 - \\beta \\frac{I_{obs}S_{obs}}{N}(1+a)(1+b) \\\\
39 \\Longrightarrow dS_{obs} &=& - \\beta \\frac{I_{obs}S_{obs}}{N}(1+b)
40 \\end{array}
42 And also:
44 .. math::
46 \\begin{array}{rcl}
47 dD &=& dD_{obs} = \\nu I = \\nu I_{obs} (1+b)
48 \\end{array}
50 And as well:
52 .. math::
54 \\begin{array}{rcl}
55 dR &=& dR_{obs} (1 + c) = \\mu I = \\mu (1 + b) I_{obs} \\\\
56 \\Longrightarrow dR_{obs} &=& - \\nu I_{obs} \\frac{1+b}{1+c}
57 \\end{array}
59 And finally:
61 .. math::
63 \\begin{array}{rcl}
64 dI &=& dI_{obs} (1 + b) = -dR - dS - dD =
65 - \\mu \\frac{1 + b}{1+ c} I_{obs} - \\nu (1+b) I_{obs} -
66 - \\beta I_{obs}\\frac{S_{obs}}{N} (1 + a)(1 + b)
67 \\\\
68 \\Longrightarrow dI_{obs} &=& - \\nu I_{obs} - \\mu I_{obs}
69 - \\beta I_{obs}\\frac{S_{obs}}{N} (1 + a)
70 \\end{array}
72 This model should still verify:
74 .. math::
76 \\begin{array}{rcl}
77 S_{obs} + I_{obs} + R_{obs} + D_{obs} &=& N = S + I + R + D \\\\
78 &=& S_{obs}(1+a) + I_{obs}(1+b) + R_{obs}(1+c) + D_{obs}
79 \\end{array}
81 That gives :math:`a S_{obs} + b I_{obs} + c R_{obs} = 0`.
83 .. runpython::
84 :showcode:
85 :rst:
87 from aftercovid.models import CovidSIRDc
89 model = CovidSIRDc()
90 print(model.to_rst())
92 .. exref::
93 :title: SIRDc simulation and plotting
95 .. plot::
97 from pandas import DataFrame
98 import matplotlib.pyplot as plt
99 from aftercovid.models import CovidSIRDc
101 model = CovidSIRDc()
102 sims = list(model.iterate(60))
103 df = DataFrame(sims)
104 print(df.head())
105 ax = df.plot(y=['S', 'I', 'R', 'D'], kind='line')
106 ax.set_xlabel("jours")
107 ax.set_ylabel("population")
108 r0 = model.R0()
109 ax.set_title("Simulation SIRDc\\nR0=%f" % r0)
111 plt.show()
113 Visual representation:
115 .. gdot::
116 :script:
118 from aftercovid.models import CovidSIRDc
119 model = CovidSIRDc()
120 print(model.to_dot())
122 See :ref:`l-base-model-sir` to get the methods
123 common to SIRx models. This model is not really working
124 better than :class:`CovidSIRD <aftercovid.covid_sird.CovidSIRD>`.
125 """
127 P0 = [
128 ('beta', 0.5, 'taux de transmission dans la population'),
129 ('mu', 1 / 14., '1/. : durée moyenne jusque la guérison'),
130 ('nu', 1 / 21., '1/. : durée moyenne jusqu\'au décès'),
131 ('a', -1.5025135094805093e-08,
132 'paramètre gérant les informations cachées (S)'),
133 ('b', 1e-5, 'paramètre gérant les informations cachées (I)'),
134 ('c', 1e-5, 'paramètre gérant les informations cachées (R)'),
135 ]
137 Q0 = [
138 ('S', 9990., 'personnes non contaminées'),
139 ('I', 10., 'nombre de personnes malades ou contaminantes'),
140 ('R', 0., 'personnes guéries (recovered)'),
141 ('D', 0., 'personnes décédées'),
142 ]
144 C0 = [
145 ('N', 10000., 'population'),
146 ]
148 # c = b (nu + mu) / (mu - b * nu)
149 # (1 + b) / (1 + c) = 1 - nu * b / mu
150 eq = {
151 'S': '-beta * (1 + b) * I * S / N',
152 'I': ('beta * (1 + a) * I * S / N'
153 '- nu * I - mu * I'),
154 'R': 'mu * (1 + b) * I / (1 + c)',
155 'D': 'nu * (1 + b) * I'}
157 def __init__(self):
158 BaseSIR.__init__(
159 self,
160 p=CovidSIRDc.P0.copy(),
161 q=CovidSIRDc.Q0.copy(),
162 c=CovidSIRDc.C0.copy(),
163 eq=CovidSIRDc.eq.copy())
165 def R0(self, t=0):
166 '''Returns R0 coefficient.
167 See :meth:`CovidSIRD.R0 <aftercovid.models.CovidSIRD.R0>`'''
168 return self['beta'] / (self['nu'] + self['mu'])
170 def correctness(self, X=None):
171 '''
172 Returns :math:`a S_{obs} + b I_{obs} + c R_{obs} = 0`.
174 :param X: None to use inner quantities
175 :return: a number
176 '''
177 if X is None:
178 X = self.vect().reshape((1, -1))
179 return (X[:, 0] * self['a'] + X[:, 1] * self['b'] +
180 X[:, 2] * self['c']) / self['N']
182 def update_abc(self, X=None, update=True, alpha=1.0, l1_ratio=0.5):
183 '''
184 Updates coefficients *a*, *b*, *c* so that method
185 :meth:`correctness <aftercovid.models.CovidSIRDc.correctness>`
186 returns 0. It uses `ElasticNet
187 <https://scikit-learn.org/stable/modules/generated/
188 sklearn.linear_model.ElasticNet.html>`_.
190 :param X: None to use inner quantities
191 :param update: True to update to the coefficients
192 or False to just return the results
193 :param alpha: see ElasticNet
194 :param l1_ratio: see ElasticNet
195 :return: dictionary
196 '''
197 if X is None:
198 X = self.vect().reshape((1, -1))
199 X = X / self['N']
200 cst = - self.correctness(X)
201 model = ElasticNet(alpha=alpha, l1_ratio=l1_ratio, fit_intercept=False)
202 model.fit(X, cst)
203 coef = model.coef_
204 res = dict(a=self['a'] + coef[0], b=self['b'] + coef[1],
205 c=self['c'] + coef[2])
206 if update:
207 self.update(**res)
208 return res
210 def rnd(self):
211 '''
212 Draws random parameters.
213 Not perfect.
214 '''
215 self['beta'] = numpy.random.randn(1) * 0.1 + 0.5
216 self['mu'] = numpy.random.randn(1) * 0.1 + 1. / 14
217 self['nu'] = numpy.random.randn(1) * 0.1 + 1. / 21
218 self['a'] = numpy.random.rand() * 1e-5
219 self['b'] = numpy.random.rand() * 1e-5
220 self['c'] = numpy.random.rand() * 1e-5
222 @staticmethod
223 def add_noise(X, epsilon=1.):
224 """
225 Tries to add reasonable noise to the quantities stored in *X*.
227 :param epsilon: amplitude
228 :return: new X
229 """
230 return CovidSIRD.add_noise(X, epsilon=epsilon)