448 views

Skip to first unread message

Oct 12, 2010, 7:45:46 PM10/12/10

to PyMC

Dear PyMCers,

I wanted to see if people had more examples of dynamic compartmental

models, such as difference equation models, monte carlo simulations

using tau-leap, etc. to estimate transition parameters.

I'm new to PyMC and was trying to do a simple SIR model to explore

it. SIR models (Susceptible-Infected-Recovered) are simple,

compartmental infectious disease models. They start with the

population being distributed in one of the three states (SIR), and

then susceptible individuals can move to the infected state at a rate

(beta) times the product of the population of infected and susceptible

individuals. Individuals then recover at a rate gamma. There are a

lot of great examples along with Python code here:

http://www.warwick.ac.uk/~masfz/ModelingInfectiousDiseases/Chapter2/Program_2.1/index.html

I want to take observed data about the population in those three

states over time to estimate the transmission parameter (beta) and

recovery parameter (gamma)

The simple SIR model can be

S[t+1] = S[t] - beta*I[t]*S[t]/N[t]

I[t+1]= I[t] + beta*I[t]*S[t]/N[t] - gamma*I[t]

R[t+1] = R[t] + gamma*I[t]

where N[t]=S[t]+I[t]+R[t]

(there are several ways to write the SIR model depending on

assumptions about transmission).

I'm a beginner to this and haven't had luck implementing it. Below is

my attempt at an SI model (without R) with made-up data and

parameters. Ideally would incorporate stochasticity at each time step,

below the 'actual data' (SIR populations) was deterministic for

simplicity.

I would be grateful for any examples, comments, or suggestions for

where to learn more. I've looked at the Fisheries Surplus model on

the PyMC site.

thank you,

Jason

##SI model

from pymc import*

import numpy as np

from numpy import *

#observed data

susceptible_array= np.array([999,997,996,994,993,992,990,989,986,984])

infected_array= np.array([1,2,5,6,7,8,9,11,13,15])

S = np.zeros(len(susceptible_array))

I = np.zeros(len(infected_array))

S[0] = 999

I[0] = 1

beta = Normal('beta', mu=0, tau=5)

gamma = Normal('gamma', mu=0, tau=5)

for i in range(1,len(susceptible_array)):

S[i] = S[i-1] - 0.05*beta*S[i-1]*I[i-1]/(S[i-1]+I[i-1])

I[i] = I[i-1] + 0.05*beta*S[i-1]*I[i-1]/(S[i-1]+I[i-1]) -

gamma*I[i-1]

A = Poisson('A', mu= S, value=susceptible_array, observed=True)

B = Poisson('B', mu= I, value=infected_array, observed=True)

I wanted to see if people had more examples of dynamic compartmental

models, such as difference equation models, monte carlo simulations

using tau-leap, etc. to estimate transition parameters.

I'm new to PyMC and was trying to do a simple SIR model to explore

it. SIR models (Susceptible-Infected-Recovered) are simple,

compartmental infectious disease models. They start with the

population being distributed in one of the three states (SIR), and

then susceptible individuals can move to the infected state at a rate

(beta) times the product of the population of infected and susceptible

individuals. Individuals then recover at a rate gamma. There are a

lot of great examples along with Python code here:

http://www.warwick.ac.uk/~masfz/ModelingInfectiousDiseases/Chapter2/Program_2.1/index.html

I want to take observed data about the population in those three

states over time to estimate the transmission parameter (beta) and

recovery parameter (gamma)

The simple SIR model can be

S[t+1] = S[t] - beta*I[t]*S[t]/N[t]

I[t+1]= I[t] + beta*I[t]*S[t]/N[t] - gamma*I[t]

R[t+1] = R[t] + gamma*I[t]

where N[t]=S[t]+I[t]+R[t]

(there are several ways to write the SIR model depending on

assumptions about transmission).

I'm a beginner to this and haven't had luck implementing it. Below is

my attempt at an SI model (without R) with made-up data and

parameters. Ideally would incorporate stochasticity at each time step,

below the 'actual data' (SIR populations) was deterministic for

simplicity.

I would be grateful for any examples, comments, or suggestions for

where to learn more. I've looked at the Fisheries Surplus model on

the PyMC site.

thank you,

Jason

##SI model

from pymc import*

import numpy as np

from numpy import *

#observed data

susceptible_array= np.array([999,997,996,994,993,992,990,989,986,984])

infected_array= np.array([1,2,5,6,7,8,9,11,13,15])

S = np.zeros(len(susceptible_array))

I = np.zeros(len(infected_array))

S[0] = 999

I[0] = 1

beta = Normal('beta', mu=0, tau=5)

gamma = Normal('gamma', mu=0, tau=5)

for i in range(1,len(susceptible_array)):

S[i] = S[i-1] - 0.05*beta*S[i-1]*I[i-1]/(S[i-1]+I[i-1])

I[i] = I[i-1] + 0.05*beta*S[i-1]*I[i-1]/(S[i-1]+I[i-1]) -

gamma*I[i-1]

A = Poisson('A', mu= S, value=susceptible_array, observed=True)

B = Poisson('B', mu= I, value=infected_array, observed=True)

Oct 14, 2010, 4:54:17 AM10/14/10

to py...@googlegroups.com

Hi Jason,

What went wrong with your attempt?

Anand

Oct 14, 2010, 8:20:06 AM10/14/10

to PyMC

On Oct 12, 6:45 pm, jason andrews <jasona...@gmail.com> wrote:

> Dear PyMCers,

>

> I wanted to see if people had more examples of dynamic compartmental

> models, such as difference equation models, monte carlo simulations

> using tau-leap, etc. to estimate transition parameters.

>

> I'm new to PyMC and was trying to do a simple SIR model to explore

> it. SIR models (Susceptible-Infected-Recovered) are simple,

> compartmental infectious disease models. They start with the

> population being distributed in one of the three states (SIR), and

> then susceptible individuals can move to the infected state at a rate

> (beta) times the product of the population of infected and susceptible

> individuals. Individuals then recover at a rate gamma. There are a

Hi Jason,

I fit some SIR models using PyMC earlier this year. They follow the

methods described in Zipkin et al:

Zipkin, E.F., Jennelle, C.S. and Cooch, E.G. 2010. A primer on the

application of Markov chains to the study of wildlife disease

dynamics. Methods in Ecology and Evolution. 1: 192-198

There are several different versions, but I will post one of them

among the examples on the PyMC project page. Look for it today.

Hope that helps,

cf

Oct 14, 2010, 2:21:11 PM10/14/10

to py...@googlegroups.com

It seems like many people had the same amount of delay in their responses. I also fit compartmental models with PyMC, and I have often wondered if this is more common than I know of, and if not, why not?

__ __

The issue with your code is you need to make the compartments a deterministic variable. I’ve cut-and-pasted your code into a github gist so that it has nice formatting: http://gist.github.com/626689

__ __

--Abie

--

You received this message because you are subscribed to the Google Groups "PyMC" group.

To post to this group, send email to py...@googlegroups.com.

To unsubscribe from this group, send email to pymc+uns...@googlegroups.com.

For more options, visit this group at http://groups.google.com/group/pymc?hl=en.

Oct 14, 2010, 7:03:42 PM10/14/10

to PyMC

Dear All,

Thanks for your quick responses and helpful comments.

Anand: When I ran the model with that code, the resultant parameters

(beta and gamma) were just giving me the distribution that I specified

rather than converging to fit the model. I didn't realize I had to

specify the compartments as deterministic parameters.

Chris: I'd love to see some more examples, that would be very

helpful. Thanks for the reference as well.

Abie: Many thanks for correcting my problem. This works perfectly. I

really appreciate your time and help.

thanks again,

Jason

On Oct 14, 2:21 pm, "A. Flaxman" <a...@uw.edu> wrote:

> It seems like many people had the same amount of delay in their responses. I also fit compartmental models with PyMC, and I have often wondered if this is more common than I know of, and if not, why not?

>

> The issue with your code is you need to make the compartments a deterministic variable. I've cut-and-pasted your code into a github gist so that it has nice formatting: http://gist.github.com/626689

>

> --Abie

>

> From: py...@googlegroups.com [mailto:py...@googlegroups.com] On Behalf Of Anand Patil

> Sent: Thursday, October 14, 2010 1:54 AM

> To: py...@googlegroups.com

> Subject: Re: [pymc] dynamic compartmental models - ex. SIR models

>

> To post to this group, send email to py...@googlegroups.com<mailto:py...@googlegroups.com>.

> To unsubscribe from this group, send email to pymc+uns...@googlegroups.com<mailto:pymc+uns...@googlegroups.com> .

Thanks for your quick responses and helpful comments.

Anand: When I ran the model with that code, the resultant parameters

(beta and gamma) were just giving me the distribution that I specified

rather than converging to fit the model. I didn't realize I had to

specify the compartments as deterministic parameters.

Chris: I'd love to see some more examples, that would be very

helpful. Thanks for the reference as well.

Abie: Many thanks for correcting my problem. This works perfectly. I

really appreciate your time and help.

thanks again,

Jason

On Oct 14, 2:21 pm, "A. Flaxman" <a...@uw.edu> wrote:

> It seems like many people had the same amount of delay in their responses. I also fit compartmental models with PyMC, and I have often wondered if this is more common than I know of, and if not, why not?

>

> The issue with your code is you need to make the compartments a deterministic variable. I've cut-and-pasted your code into a github gist so that it has nice formatting: http://gist.github.com/626689

>

> --Abie

>

> From: py...@googlegroups.com [mailto:py...@googlegroups.com] On Behalf Of Anand Patil

> Sent: Thursday, October 14, 2010 1:54 AM

> To: py...@googlegroups.com

> Subject: Re: [pymc] dynamic compartmental models - ex. SIR models

>

> On Wed, Oct 13, 2010 at 12:45 AM, jason andrews <jasona...@gmail.com<mailto:jasona...@gmail.com>> wrote:

> Dear PyMCers,

>

> I wanted to see if people had more examples of dynamic compartmental

> models, such as difference equation models, monte carlo simulations

> using tau-leap, etc. to estimate transition parameters.

>

> I'm new to PyMC and was trying to do a simple SIR model to explore

> it. SIR models (Susceptible-Infected-Recovered) are simple,

> compartmental infectious disease models. They start with the

> population being distributed in one of the three states (SIR), and

> then susceptible individuals can move to the infected state at a rate

> (beta) times the product of the population of infected and susceptible

> individuals. Individuals then recover at a rate gamma. There are a

> lot of great examples along with Python code here:http://www.warwick.ac.uk/~masfz/ModelingInfectiousDiseases/Chapter2/P...
> Dear PyMCers,

>

> I wanted to see if people had more examples of dynamic compartmental

> models, such as difference equation models, monte carlo simulations

> using tau-leap, etc. to estimate transition parameters.

>

> I'm new to PyMC and was trying to do a simple SIR model to explore

> it. SIR models (Susceptible-Infected-Recovered) are simple,

> compartmental infectious disease models. They start with the

> population being distributed in one of the three states (SIR), and

> then susceptible individuals can move to the infected state at a rate

> (beta) times the product of the population of infected and susceptible

> individuals. Individuals then recover at a rate gamma. There are a

> To unsubscribe from this group, send email to pymc+uns...@googlegroups.com<mailto:pymc+uns...@googlegroups.com> .

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu