search

mrc-ide / odin

ᚩ A DSL for describing and solving differential equations in R
https://mrc-ide.github.io/odin
Other
106 stars 13 forks source link

How to replace SIR Model and use PMCMC to estimate parameters #307

Closed whzhuscu closed 1 year ago

whzhuscu commented 1 year ago

What is code of the example "sir" model?

gen_sir <- dust::dust_example("sir")

When I use the code to write sir model:

library(odin.dust)
gen_sir <- odin.dust::odin_dust("sir.R")
ℹ 20 functions decorated with [[cpp11::register]]
✔ generated file 'cpp11.R'
✔ generated file 'cpp11.cpp'
Re-compiling sir76688fd8

─  installing *source* package 'sir76688fd8' ... (570ms)

   ** using staged installation
   ** libs

   "D:/Program Files/R/rtools/rtools40/mingw64/bin/"g++ -std=gnu++11  -I"D:/PROGRA~1/R/R-42~1.1/include" -DNDEBUG  -I'D:/Program Files/R/R-4.2.1/library/cpp11/include'     -IC:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include -DHAVE_INLINE -fopenmp    -O2 -Wall  -mfpmath=sse -msse2 -mstackrealign -c cpp11.cpp -o cpp11.o

   "D:/Program Files/R/rtools/rtools40/mingw64/bin/"g++ -std=gnu++11  -I"D:/PROGRA~1/R/R-42~1.1/include" -DNDEBUG  -I'D:/Program Files/R/R-4.2.1/library/cpp11/include'     -IC:/Users/Administrator/AppData/Local/R/win-library/4.2/dust/include -DHAVE_INLINE -fopenmp    -O2 -Wall  -mfpmath=sse -msse2 -mstackrealign -c dust.cpp -o dust.o

   D:/Program Files/R/rtools/rtools40/mingw64/bin/g++ -std=gnu++11 -shared -s -static-libgcc -o sir76688fd8.dll tmp.def cpp11.o dust.o -fopenmp -LD:/PROGRA~1/R/R-42~1.1/bin/x64 -lR

   installing to C:/Users/ADMINI~1/AppData/Local/Temp/Rtmpey6Alb/devtools_install_4cfc1a634851/00LOCK-file4cfc563a3e37/00new/sir76688fd8/libs/x64

─  DONE (sir76688fd8)

ℹ Loading sir76688fd8

there is error means "subscript out of bounds":

> incidence <- read.csv(system.file("sir_incidence.csv", package = "mcstate"))
> dt <- 0.25
> sir_data <- mcstate::particle_filter_data(data = incidence,
+                                           time = "day",
+                                           rate = 1 / dt)
Warning in mcstate::particle_filter_data(data = incidence, time = "day",  :
  'initial_time' should be provided. I'm assuming '0' which is one time unit before the first time in your data (1), but this might not be appropriate. This will become an error in a future version of mcstate
> rmarkdown::paged_table(sir_data)
    day_start day_end time_start time_end cases
1           0       1          0        4     3
2           1       2          4        8     2
3           2       3          8       12     2
4           3       4         12       16     2
5           4       5         16       20     1
6           5       6         20       24     3
7           6       7         24       28     2
8           7       8         28       32     5
9           8       9         32       36     5
10          9      10         36       40     6
11         10      11         40       44     9
12         11      12         44       48     3
13         12      13         48       52     5
14         13      14         52       56     7
15         14      15         56       60     9
16         15      16         60       64    12
17         16      17         64       68    13
18         17      18         68       72    13
19         18      19         72       76    17
20         19      20         76       80    16
21         20      21         80       84    14
22         21      22         84       88    24
23         22      23         88       92    16
24         23      24         92       96    16
25         24      25         96      100    21
26         25      26        100      104    18
27         26      27        104      108    17
28         27      28        108      112    14
29         28      29        112      116    15
30         29      30        116      120    18
31         30      31        120      124    14
32         31      32        124      128    13
33         32      33        128      132    17
34         33      34        132      136     9
35         34      35        136      140    21
36         35      36        140      144    15
37         36      37        144      148    20
38         37      38        148      152    19
39         38      39        152      156    18
40         39      40        156      160    20
41         40      41        160      164    21
42         41      42        164      168    25
43         42      43        168      172    14
44         43      44        172      176     7
45         44      45        176      180    13
46         45      46        180      184     9
47         46      47        184      188    18
48         47      48        188      192    12
49         48      49        192      196    10
50         49      50        196      200    16
51         50      51        200      204     6
52         51      52        204      208    10
53         52      53        208      212     4
54         53      54        212      216    15
55         54      55        216      220     3
56         55      56        220      224     4
57         56      57        224      228     8
58         57      58        228      232     8
59         58      59        232      236     2
60         59      60        236      240     2
61         60      61        240      244     5
62         61      62        244      248     5
63         62      63        248      252     4
64         63      64        252      256     3
65         64      65        256      260     2
66         65      66        260      264     5
67         66      67        264      268     6
68         67      68        268      272     0
69         68      69        272      276     3
70         69      70        276      280     3
71         70      71        280      284     4
72         71      72        284      288     1
73         72      73        288      292     1
74         73      74        292      296     0
75         74      75        296      300     3
76         75      76        300      304     3
77         76      77        304      308     4
78         77      78        308      312     1
79         78      79        312      316     3
80         79      80        316      320     1
81         80      81        320      324     0
82         81      82        324      328     0
83         82      83        328      332     3
84         83      84        332      336     3
85         84      85        336      340     1
86         85      86        340      344     0
87         86      87        344      348     0
88         87      88        348      352     3
89         88      89        352      356     1
90         89      90        356      360     3
91         90      91        360      364     7
92         91      92        364      368     5
93         92      93        368      372     0
94         93      94        372      376     1
95         94      95        376      380     0
96         95      96        380      384     4
97         96      97        384      388     1
98         97      98        388      392     0
99         98      99        392      396     5
100        99     100        396      400     1
> plot(incidence$day, incidence$cases,
+      type = "l", xlab = "Day", ylab = "New cases")

> #### Defining the comparison function ####
> case_compare <- function(state, observed, pars = NULL) {
+   exp_noise <- 1e6
+   
+   incidence_modelled <- state[5, , drop = TRUE]
+   incidence_observed <- observed$cases
+   lambda <- incidence_modelled +
+     rexp(n = length(incidence_modelled), rate = exp_noise)
+   dpois(x = incidence_observed, lambda = lambda, log = TRUE)
+ }
> gen_sir$new(pars = list(), time = 0, n_particles = 1L)$info()
$dim
$dim$time
[1] 1

$dim$S
[1] 1

$dim$I
[1] 1

$dim$R
[1] 1

$len
[1] 4

$index
$index$time
[1] 1

$index$S
[1] 2

$index$I
[1] 3

$index$R
[1] 4
> incidence_compare <- function(state, prev_state, observed, pars = NULL) {
+   exp_noise <- 1e6
+   
+   lambda <- state[4, , drop = TRUE] +
+     rexp(n = length(incidence_modelled), rate = exp_noise)
+   dpois(x = observed$cases, lambda = lambda, log = TRUE)
+ }

> #### Inferring parameters ####
> n_particles <- 100
> filter <- mcstate::particle_filter$new(data = sir_data,
+                                        model = gen_sir,
+                                        n_particles = n_particles,
+                                        compare = case_compare,
+                                        seed = 1L)
> filter$run(save_history = TRUE, pars = list(dt = dt))
Error in state[5, , drop = TRUE] : Subscript out of bounds
richfitz commented 1 year ago

The answer remains the same as the last time you asked this here: https://github.com/mrc-ide/odin/issues/304

Please see the dust documentation, where this is discussed: mrc-ide.github.io/dust/articles/dust.html

More documentation around the set of packages can be found here: mrc-ide.github.io/odin/articles/guide.html

whzhuscu commented 1 year ago

Thanks!