inst/examples/gillespie.sir.R

In pomp: Statistical Inference for Partially Observed Markov Processes

require(pomp)

dat <- '"time";"reports"
0;0
0,0192307692307692;5
0,0384615384615385;5
0,0576923076923077;2
0,0769230769230769;2
0,0961538461538462;3
0,115384615384615;2
0,134615384615385;2
0,153846153846154;1
0,173076923076923;3
0,192307692307692;2
0,211538461538462;3
0,230769230769231;1
0,25;1
0,269230769230769;2
0,288461538461538;1
0,307692307692308;2
0,326923076923077;3
0,346153846153846;1
0,365384615384615;0
0,384615384615385;0
0,403846153846154;2
0,423076923076923;1
0,442307692307692;0
0,461538461538462;1
0,480769230769231;2
0,5;0
0,519230769230769;0
0,538461538461539;0
0,557692307692308;1
0,576923076923077;2
0,596153846153846;0
0,615384615384615;1
0,634615384615385;2
0,653846153846154;0
0,673076923076923;2
0,692307692307692;2
0,711538461538462;3
0,730769230769231;1
0,75;1
0,769230769230769;0
0,788461538461539;1
0,807692307692308;1
0,826923076923077;0
0,846153846153846;3
0,865384615384615;2
0,884615384615385;10
0,903846153846154;3
0,923076923076923;4
0,942307692307692;2
0,961538461538462;2
0,980769230769231;11
1;4
1,01923076923077;8
1,03846153846154;4
1,05769230769231;1
1,07692307692308;7
1,09615384615385;7
1,11538461538462;5
1,13461538461538;3
1,15384615384615;2
1,17307692307692;7
1,19230769230769;9
1,21153846153846;9
1,23076923076923;5
1,25;3
1,26923076923077;4
1,28846153846154;12
1,30769230769231;3
1,32692307692308;11
1,34615384615385;2
1,36538461538462;10
1,38461538461538;9
1,40384615384615;5
1,42307692307692;10
1,44230769230769;10
1,46153846153846;8
1,48076923076923;7
1,5;7
1,51923076923077;5
1,53846153846154;10
1,55769230769231;10
1,57692307692308;13
1,59615384615385;14
1,61538461538462;22
1,63461538461538;18
1,65384615384615;22
1,67307692307692;16
1,69230769230769;27
1,71153846153846;24
1,73076923076923;24
1,75;36
1,76923076923077;42
1,78846153846154;35
1,80769230769231;43
1,82692307692308;50
1,84615384615385;62
1,86538461538462;47
1,88461538461538;58
1,90384615384615;50
1,92307692307692;57
1,94230769230769;51
1,96153846153846;50
1,98076923076923;54
2;61
2,01923076923077;64
2,03846153846154;58
2,05769230769231;61
2,07692307692308;80
2,09615384615385;48
2,11538461538462;61
2,13461538461538;59
2,15384615384615;59
2,17307692307692;57
2,19230769230769;39
2,21153846153846;42
2,23076923076923;45
2,25;55
2,26923076923077;26
2,28846153846154;40
2,30769230769231;50
2,32692307692308;40
2,34615384615385;49
2,36538461538462;44
2,38461538461538;57
2,40384615384615;44
2,42307692307692;56
2,44230769230769;40
2,46153846153846;54
2,48076923076923;59
2,5;67
2,51923076923077;52
2,53846153846154;51
2,55769230769231;57
2,57692307692308;66
2,59615384615385;69
2,61538461538462;61
2,63461538461538;63
2,65384615384615;55
2,67307692307692;59
2,69230769230769;59
2,71153846153846;75
2,73076923076923;62
2,75;66
2,76923076923077;80
2,78846153846154;53
2,80769230769231;69
2,82692307692308;53
2,84615384615385;56
2,86538461538462;50
2,88461538461538;45
2,90384615384615;41
2,92307692307692;40
2,94230769230769;32
2,96153846153846;39
2,98076923076923;35
3;25
3,01923076923077;32
3,03846153846154;15
3,05769230769231;20
3,07692307692308;19
3,09615384615385;16
3,11538461538462;11
3,13461538461538;14
3,15384615384615;16
3,17307692307692;11
3,19230769230769;8
3,21153846153846;6
3,23076923076923;3
3,25;8
3,26923076923077;8
3,28846153846154;2
3,30769230769231;3
3,32692307692308;4
3,34615384615385;8
3,36538461538462;7
3,38461538461539;4
3,40384615384615;9
3,42307692307692;7
3,44230769230769;8
3,46153846153846;5
3,48076923076923;5
3,5;7
3,51923076923077;1
3,53846153846154;6
3,55769230769231;7
3,57692307692308;6
3,59615384615385;5
3,61538461538462;8
3,63461538461539;10
3,65384615384615;1
3,67307692307692;0
3,69230769230769;5
3,71153846153846;8
3,73076923076923;1
3,75;8
3,76923076923077;8
3,78846153846154;5
3,80769230769231;7
3,82692307692308;3
3,84615384615385;7
3,86538461538462;4
3,88461538461539;6
3,90384615384615;7
3,92307692307692;6
3,94230769230769;7
3,96153846153846;7
3,98076923076923;5
4;5
4,01923076923077;7
4,03846153846154;2
4,05769230769231;8
4,07692307692308;3
4,09615384615385;3
4,11538461538462;5
4,13461538461539;1
4,15384615384615;3
4,17307692307692;5
4,19230769230769;5
4,21153846153846;2
4,23076923076923;3
4,25;6
4,26923076923077;1
4,28846153846154;2
4,30769230769231;2
4,32692307692308;5
4,34615384615385;3
4,36538461538462;3
4,38461538461539;7
4,40384615384615;2
4,42307692307692;2
4,44230769230769;2
4,46153846153846;3
4,48076923076923;4
4,5;4
4,51923076923077;7
4,53846153846154;1
4,55769230769231;5
4,57692307692308;12
4,59615384615385;5
4,61538461538462;8
4,63461538461539;3
4,65384615384615;17
4,67307692307692;14
4,69230769230769;15
4,71153846153846;12
4,73076923076923;19
4,75;19
4,76923076923077;17
4,78846153846154;23
4,80769230769231;28
4,82692307692308;23
4,84615384615385;32
4,86538461538462;24
4,88461538461539;35
4,90384615384615;31
4,92307692307692;23
4,94230769230769;41
4,96153846153846;27
4,98076923076923;28
5;26
5,01923076923077;29
5,03846153846154;38
5,05769230769231;20
5,07692307692308;17
5,09615384615385;29
5,11538461538462;31
5,13461538461539;25
5,15384615384615;31
5,17307692307692;26
5,19230769230769;23
5,21153846153846;36
5,23076923076923;25
5,25;25
5,26923076923077;35
5,28846153846154;37
5,30769230769231;32
5,32692307692308;30
5,34615384615385;41
5,36538461538462;45
5,38461538461539;61
5,40384615384615;40
5,42307692307692;43
5,44230769230769;46
5,46153846153846;49
5,48076923076923;59
5,5;48
5,51923076923077;69
5,53846153846154;61
5,55769230769231;68
5,57692307692308;79
5,59615384615385;91
5,61538461538462;93
5,63461538461539;77
5,65384615384615;87
5,67307692307692;107
5,69230769230769;103
5,71153846153846;119
5,73076923076923;111
5,75;110
5,76923076923077;91
5,78846153846154;120
5,80769230769231;106
5,82692307692308;93
5,84615384615385;69
5,86538461538462;90
5,88461538461539;100
5,90384615384615;73
5,92307692307692;72
5,94230769230769;67
5,96153846153846;63
5,98076923076923;41
6;43
6,01923076923077;33
6,03846153846154;40
6,05769230769231;29
6,07692307692308;25
6,09615384615385;22
6,11538461538462;19
6,13461538461539;15
6,15384615384615;18
6,17307692307692;13
6,19230769230769;15
6,21153846153846;14
6,23076923076923;10
6,25;8
6,26923076923077;13
6,28846153846154;4
6,30769230769231;10
6,32692307692308;3
6,34615384615385;9
6,36538461538462;10
6,38461538461539;9
6,40384615384615;4
6,42307692307692;4
6,44230769230769;3
6,46153846153846;2
6,48076923076923;6
6,5;4
6,51923076923077;5
6,53846153846154;3
6,55769230769231;8
6,57692307692308;3
6,59615384615385;6
6,61538461538462;5
6,63461538461539;6
6,65384615384615;1
6,67307692307692;8
6,69230769230769;12
6,71153846153846;10
6,73076923076923;3
6,75;3
6,76923076923077;6
6,78846153846154;5
6,80769230769231;10
6,82692307692308;6
6,84615384615385;6
6,86538461538462;5
6,88461538461539;11
6,90384615384615;6
6,92307692307692;14
6,94230769230769;7
6,96153846153846;6
6,98076923076923;5
7;4
7,01923076923077;7
7,03846153846154;3
7,05769230769231;4
7,07692307692308;3
7,09615384615385;3
7,11538461538462;2
7,13461538461539;3
7,15384615384615;7
7,17307692307692;2
7,19230769230769;2
7,21153846153846;3
7,23076923076923;3
7,25;5
7,26923076923077;2
7,28846153846154;3
7,30769230769231;7
7,32692307692308;2
7,34615384615385;6
7,36538461538462;2
7,38461538461539;5
7,40384615384615;5
7,42307692307692;6
7,44230769230769;5
7,46153846153846;8
7,48076923076923;9
7,5;5
7,51923076923077;5
7,53846153846154;15
7,55769230769231;9
7,57692307692308;18
7,59615384615385;17
7,61538461538462;10
7,63461538461539;16
7,65384615384615;15
7,67307692307692;19
7,69230769230769;14
7,71153846153846;23
7,73076923076923;28
7,75;32
7,76923076923077;32
7,78846153846154;31
7,80769230769231;40
7,82692307692308;37
7,84615384615385;41
7,86538461538462;37
7,88461538461539;40
7,90384615384615;38
7,92307692307692;44
7,94230769230769;35
7,96153846153846;34
7,98076923076923;26
8;49
8,01923076923077;32
8,03846153846154;45
8,05769230769231;38
8,07692307692308;31
8,09615384615385;26
8,11538461538461;27
8,13461538461539;26
8,15384615384616;26
8,17307692307692;22
8,19230769230769;33
8,21153846153846;36
8,23076923076923;27
8,25;36
8,26923076923077;39
8,28846153846154;36
8,30769230769231;30
8,32692307692308;39
8,34615384615385;28
8,36538461538461;31
8,38461538461539;40
8,40384615384616;47
8,42307692307692;41
8,44230769230769;42
8,46153846153846;42
8,48076923076923;50
8,5;60
8,51923076923077;53
8,53846153846154;64
8,55769230769231;59
8,57692307692308;65
8,59615384615385;61
8,61538461538462;69
8,63461538461539;66
8,65384615384616;71
8,67307692307692;71
8,69230769230769;72
8,71153846153846;94
8,73076923076923;79
8,75;79
8,76923076923077;78
8,78846153846154;78
8,80769230769231;92
8,82692307692308;64
8,84615384615385;70
8,86538461538462;50
8,88461538461539;63
8,90384615384616;58
8,92307692307692;53
8,94230769230769;48
8,96153846153846;51
8,98076923076923;28
9;36
9,01923076923077;40
9,03846153846154;20
9,05769230769231;16
9,07692307692308;19
9,09615384615385;16
9,11538461538462;20
9,13461538461539;16
9,15384615384616;11
9,17307692307692;9
9,19230769230769;15
9,21153846153846;15
9,23076923076923;7
9,25;12
9,26923076923077;9
9,28846153846154;5
9,30769230769231;4
9,32692307692308;3
9,34615384615385;4
9,36538461538462;4
9,38461538461539;5
9,40384615384616;3
9,42307692307692;2
9,44230769230769;2
9,46153846153846;3
9,48076923076923;3
9,5;5
9,51923076923077;3
9,53846153846154;3
9,55769230769231;6
9,57692307692308;6
9,59615384615385;5
9,61538461538462;2
9,63461538461539;5
9,65384615384616;5
9,67307692307692;3
9,69230769230769;3
9,71153846153846;3
9,73076923076923;4
9,75;3
9,76923076923077;9
9,78846153846154;4
9,80769230769231;2
9,82692307692308;5
9,84615384615385;4
9,86538461538462;2
9,88461538461539;4
9,90384615384616;5
9,92307692307692;3
9,94230769230769;2
9,96153846153846;3
9,98076923076923;2
10;6
'

pomp(
     data=read.csv2(text=dat),
     times="time",
     t0=0,
     params=c(
       gamma=24,mu=1/70,iota=0.1,
       beta1=330,beta2=410,beta3=490,
       rho=0.1,
       S_0=0.05,I_0=1e-4,R_0=0.95,
       pop=1000000,
       beta.sd=0
       ),
     rprocess=gillespie.sim(
       rate.fun="_sir_rates",
       PACKAGE="pomp",
       v=cbind(
         birth=c(1,0,0,1,0),
         sdeath=c(-1,0,0,-1,0),
         infection=c(-1,1,0,0,0),
         ideath=c(0,-1,0,-1,0),
         recovery=c(0,-1,1,0,1),
         rdeath=c(0,0,-1,-1,0)
         ),
       d=cbind(
         birth=c(0,0,0,1,0),
         sdeath=c(1,0,0,0,0),
         infection=c(1,1,0,1,0),
         ideath=c(0,1,0,0,0),
         recovery=c(0,1,0,0,0),
         rdeath=c(0,0,1,0,0)
         )
       ),
     skeleton.type="vectorfield",
     skeleton="_sir_ODE",
     measurement.model=reports~binom(size=cases,prob=rho),
     PACKAGE="pomp",
     statenames=c("S","I","R","N","cases"),
     paramnames=c(
       "gamma","mu","iota",
       "beta1","beta.sd","pop","rho",
       "S_0","I_0","R_0"
       ),
     zeronames=c("cases"),
     comp.names=c("S","I","R"),
     ic.names=c("S_0","I_0","R_0"),
     fromEstimationScale="_sir_par_trans",
     toEstimationScale="_sir_par_untrans",
     nbasis=3L,
     degree=3L,
     period=1.0,
     initializer=function(params, t0, comp.names, ic.names, ...) {
       x0 <- setNames(numeric(5),c("S","I","R","N","cases"))
       fracs <- params[ic.names]
       x0["N"] <- params["pop"]
       x0[comp.names] <- round(params["pop"]*fracs/sum(fracs))
       x0
     }
     ) -> gillespie.sir

## originally, the data were created via:
## simulate(po,nsim=1,seed=1165270654L) -> gillespie.sir

c("gillespie.sir")