SIR model with Python

Поделиться
HTML-код
  • Опубликовано: 15 янв 2025
  • Example for SIR model with Python. (Jupyter Notebook.)
  • НаукаНаука

Комментарии • 54

  • @JeanDAVID
    @JeanDAVID 4 года назад +26

    import numpy as np
    import scipy.integrate
    import matplotlib.pyplot as plt
    def SIRD_model(y, t, beta, gamma, delta): #adding death separately from recovery
    S, I, R, D = y
    dS = -beta*S*I
    dI = beta*S*I - gamma*I - delta*I
    dR = gamma*I
    dD = delta*I
    return ([dS, dI, dR, dD])
    t = np.linspace(0, 100,1000)
    S0 = 0.9
    I0 = 0.1
    R0 = 0.0
    D0 = 0.0
    beta = 0.45
    gamma = 0.1
    delta = 0.15
    Ro = beta/(delta+gamma)
    print('R0 =',Ro)
    solution = scipy.integrate.odeint(SIRD_model, [S0, I0, R0, D0], t, args=(beta, gamma, delta))
    plt.figure(figsize=[6,4])
    plt.plot(t, solution[:, 0], label="S(t)")
    plt.plot(t, solution[:, 1], label="I(t)")
    plt.plot(t, solution[:, 2], label="R(t)")
    plt.plot(t, solution[:, 3], label="D(t)")
    plt.grid()
    plt.legend()
    plt.xlabel("Time")
    plt.ylabel("Proportions")
    plt.title("SIRD model")
    plt.show()

  • @robertmunroe9635
    @robertmunroe9635 4 года назад +3

    ah you got me. its not linespace its linspace.

  • @bartekburmistrz8679
    @bartekburmistrz8679 8 месяцев назад

    Thank you, I have this due in 18 minutes

  • @klam77
    @klam77 4 года назад

    Beautiful

  • @ezekielizedonmi310
    @ezekielizedonmi310 4 года назад +1

    Pls i do not understand the S,I,R = y. what is 'y'? ------ SIR_model function

  • @chamba149
    @chamba149 4 года назад

    Is the time measured in days?

  • @basic-games621
    @basic-games621 4 года назад

    How can I model the exact same but instead of graphs to me like point ?
    And How can I add E=exposed into this code ?

  • @robertkavensky9709
    @robertkavensky9709 7 лет назад +2

    Hello i'm looking for the computer code of SIR for python in order to do some experiments. Could you please send me your code ( program)? It will be really nice.

    • @computationalscientist6368
      @computationalscientist6368  7 лет назад +1

      Hi! I'm afraid I don't have the code anymore, but you can easily copy it from the video which can't take more than 1-2 minutes. I hope you found this tutorial helpful.

    • @miguelpuentes9643
      @miguelpuentes9643 6 лет назад

      sorry i get an error, line 29, in
      solution = scipy.integrate.odeint(modelo_sir, [S0, I0, R0], t, args=(beta, gamma))
      File "C:\Users\HP\AppData\Local\Programs\Python\Python36\lib\site-packages\scipy\integrate\odepack.py", line 233, in odeint
      int(bool(tfirst)))
      RuntimeError: The size of the array returned by func (2) does not match the size of y0 (3).

    • @computationalscientist6368
      @computationalscientist6368  6 лет назад

      Hi! I'm afraid I can't figure out the problem from the error message. Can you please check again that your code is the same as the code in the video? If it is the same, can you post your code here? I will have time to check it during the evening (EU time).

    • @miguelpuentes9643
      @miguelpuentes9643 6 лет назад +2

      import scipy.integrate
      import numpy
      import matplotlib.pyplot as plt
      def modelo_sir(y, t, beta, gamma):
      S, I, R = y
      dS_dt = -beta*S*I
      dI_dt = beta*S*I - gamma*I
      dR_dt = gamma*I
      return([dS_dt, dI_dt,])
      # Condiciones iniciales
      S0 = 0.9
      I0 = 0.1
      R0 = 0.0
      beta = 0.3
      gamma = 0.1
      # Vector tiempo
      t = numpy.linspace(0, 100, 1000)
      # Resultado
      solucion = scipy.integrate.odeint(modelo_sir, [S0, I0, R0], t, args=(beta, gamma))
      solucion = numpy.array(solucion)
      # Resultado en grafica
      plt.figure(figsize=[6, 4])
      plt.plot(t, solucion[:, 0], label="S(t)")
      plt.plot(t, solucion[:, 1], label="I(t)")
      plt.plot(t, solucion[:, 2], label="R(t)")
      plt.grid()
      plt.legend()
      plt.xlabel("Tiempo")
      plt.ylabel("Proporciones")
      plt.title("Modelo SIR (Epidemia)")
      plt.show()

    • @miguelpuentes9643
      @miguelpuentes9643 6 лет назад +1

      I already solved the problem, thanks, like

  • @AlessandroBottoni
    @AlessandroBottoni 4 года назад +6

    For the people interested in source code, here there is a implementation of the basic SIR model in Python using SciPy (NumPy and other packages): scipython.com/book/chapter-8-scipy/additional-examples/the-sir-epidemic-model/ . Here another one, from the epirecipes project: epirecip.es/epicookbook/chapters/kr08/2_1/python_original . Searching "sir model python" on the web, you can find many others examples.

    • @youupiiiyouupiii8235
      @youupiiiyouupiii8235 4 года назад

      Hello, in your publication you say that gamma represents the inverse of the average number of days to heal. In this case, if I set beta to zero, there will be no more new infected people and the number of infected people will have to go to zero after ten days plus at least a certain interval since it's an average. I've done the test, the number of infected people goes to zero around the fiftieth day. Can you please explain to me what is wrong with my thinking?

    • @marbellydavila3065
      @marbellydavila3065 4 года назад

      I tried to run the code you shared above and I am getting errors when running the last part "Plotting the data" . :-(

    • @computationalscientist6368
      @computationalscientist6368  4 года назад

      Can you share your code?

    • @México_Maravilloso
      @México_Maravilloso 4 года назад +1

      @@marbellydavila3065
      import numpy as np
      from scipy.integrate import odeint
      import matplotlib.pyplot as plt
      # Total population, N.
      N = 1000
      # 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 = 0.2, 1./10
      # A grid of time points (in days)
      t = np.linspace(0, 150, 150)
      # 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 the time grid, t.
      ret = odeint(deriv, y0, t, args=(N, beta, gamma))
      S, I, R = ret.T
      # Plot the data on three separate curves for S(t), I(t) and R(t)
      fig = plt.figure(facecolor='w')
      ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)
      ax.plot(t, S/1000, 'b', alpha=0.5, lw=2, label='Susceptible')
      ax.plot(t, I/1000, 'r', alpha=0.5, lw=2, label='Infected')
      ax.plot(t, R/1000, 'g', alpha=0.5, lw=2, label='Recovered with immunity')
      ax.set_xlabel('Time /days')
      ax.set_ylabel('Number (1000s)')
      ax.set_ylim(0,1.2)
      ax.yaxis.set_tick_params(length=0)
      ax.xaxis.set_tick_params(length=0)
      ax.grid(b=True, which='major', c='w', lw=2, ls='-')
      legend = ax.legend()
      legend.get_frame().set_alpha(0.5)
      for spine in ('top', 'right', 'bottom', 'left'):
      ax.spines[spine].set_visible(False)
      plt.show()

    • @México_Maravilloso
      @México_Maravilloso 4 года назад

      ax = fig.add_subplot(111, facecolor='#dddddd', axisbelow=True)

  • @computersense5869
    @computersense5869 3 года назад

    what is S0

  • @yassine_pastis53
    @yassine_pastis53 2 года назад

    what is the y in S,I,R=y

  • @Dafiya445
    @Dafiya445 2 года назад

    Thank You Somuch

  • @dishanaik6128
    @dishanaik6128 3 года назад

    Thank you so much...

  • @KaushalSoni2205
    @KaushalSoni2205 Год назад

    Make a vedio for bifurcation diagram

  • @kaoutarbentalha7620
    @kaoutarbentalha7620 4 года назад

    how can we apply that model on a dataset(csv file format)?

    • @computationalscientist6368
      @computationalscientist6368  4 года назад +2

      Would you like to calibrate the model to an existing dataset?

    • @pedrosierra6102
      @pedrosierra6102 4 года назад

      @@computationalscientist6368 that is exactly what i am looking for, can you make a video to fit data on this ODE system? Please

    • @klam77
      @klam77 4 года назад

      Apply regression

    • @kaoutarbentalha7620
      @kaoutarbentalha7620 4 года назад

      @@klam77 i mean to determine those.coefficient from our dataset?

    • @klam77
      @klam77 4 года назад +1

      @@kaoutarbentalha7620 maybe i'm missing something (you can tell me). (1) you have csv with infection counts. (2) you have SIR model, where I(t) = exp(aS-b)*t. (3) you must FIT your csv infection count data to the exponential: (4) take the log of Infection count and (5) run a linear regression vs time. (6) this will give you (aS-b) coefficient in toto, and (7) then use eqn for R(t) = b*exp(as-b)*t and deduce the "b". (8) then you can back out the "a" from earlier, now that you know "b";
      thusly you have fitted the entire SIR model...no?

  • @bryanjacobo4963
    @bryanjacobo4963 4 года назад

    Why R0 is 0 if in the SIR model R0 is beta/gamma?

    • @thusithawijesooriya1463
      @thusithawijesooriya1463 4 года назад

      Initial Removed value

    • @zhangsc91
      @zhangsc91 4 года назад

      Yes, unfortunate notation but initial removed R(0) vs basic reproduction R_0 are two very different things..

  • @lilymunro9914
    @lilymunro9914 4 года назад

    is there a way to change the step size of t in this model?

    • @zhangsc91
      @zhangsc91 4 года назад

      linspace: 3rd argument controls how many total time steps to take

  • @JuanFernandoAldanaWong
    @JuanFernandoAldanaWong 4 года назад +1

    ¡No sirve!