-
Notifications
You must be signed in to change notification settings - Fork 3
/
Example2.R
35 lines (27 loc) · 1.5 KB
/
Example2.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
# data from UW epidemic as described in H1N1 paper with ES
library(plyr)
library(MASS)
library(smfsb)
library(ramcmc)
library(SDSMCMC)
############################################################################
# Second example data: simulation data. The data is simulated by Sellke construction.
# This example has two MCMC simulations using synthetic data by Sellke construction.
#The first is estimation including estimation n, number of susceptible. The second runs without n.
burn = 1000; # burning period
thin =1; # tinning period
nrepeat = 2000; # number of posterior sample;
sim.num= burn + nrepeat * thin; # total number of simulation
#initial parameter setting
k1 = 1.1; k2 = 0.8 ; k3 = 0.05; n=1000; T.max = 30;
beta=k1; gamma=k2; rho=k3;
#generating synthetic epidemic data using Sellke construction
pop.data = Sellke(n=n, rho=k3, beta=k1, gamma=k2, Tmax = T.max)
#MCMC using Method 3. In this example, we estimate three parameters given n = 1000
sds <- SDS.MCMC(data = pop.data, Tmax=T.max, fitn = F, nrepeat = sim.num,
prior.a=c(0.001,0.001,0.001), prior.b=c(0.001,0.001,0.001), ic = c(k1, k2, k3, n))
result(data = sds, fitn = F)
# MCMC using Method 3. In this example, we estimate three parameters and the initial number of susecptible
sds <- SDS.MCMC(data = pop.data, Tmax=T.max, fitn = T, nrepeat = sim.num,
prior.a=c(0.001,0.001,0.001), prior.b=c(0.001,0.001,0.001), ic = c(k1, k2, k3, n))
result(data = sds, fitn = T)