Simple prediction using Katla data alone...
model {
for (i in 2:(Nk-1)) {
t[i] <- D[i] - D[i-1]
}
for (i in 2:Nk) {
t[i] ~ dweib(r, mu)I(cens[i],)
}
r ~ dunif(0, 10)
sigma ~ dgamma(0.001, 0.001) # Jeffreys prior
mu <- 1/pow(sigma, r)
median <- sigma*pow(log(2), 1/r) # median time to event
p.erupt.1 <- 1 - exp(pow(92/sigma,r) - pow((92+1)/sigma,r))
p.erupt.5 <- 1 - exp(pow(92/sigma,r) - pow((92+5)/sigma,r))
p.erupt.10 <- 1 - exp(pow(92/sigma,r) - pow((92+10)/sigma,r))
p.erupt.50 <- 1 - exp(pow(92/sigma,r) - pow((92+50)/sigma,r))
}
Data:
list(Nk=19,
D = c(1177, 1262, 1311, 1357, 1416, 1440, 1450, 1500, 1550,
1580, 1612, 1625, 1660, 1721, 1755, 1823, 1860, 1918, NA),
cens = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,92))
Inits:
list(r=1, sigma=10)
node mean sd MC error 2.5% median 97.5% start sample
t[19] 106.2 16.22 0.3188 92.33 101.3 148.4 1000 5001
median 45.23 6.158 0.1093 33.42 45.17 57.87 1000 5001
p.erupt.1 0.07364 0.03261 5.453E-4 0.02759 0.06752 0.1525 1000 5001
p.erupt.5 0.3161 0.1157 0.002027 0.132 0.3009 0.5782 1000 5001
p.erupt.10 0.5286 0.1539 0.002845 0.2502 0.5212 0.8376 1000 5001
p.erupt.50 0.9647 0.05829 0.00123 0.7927 0.9892 1.0 1000 5001
r 2.111 0.4266 0.008172 1.376 2.085 3.028 1000 5001
sigma 54.13 6.865 0.1215 41.9 53.85 68.91 1000 5001
mu 7.828E-4 0.00176 3.57E-5 4.612E-6 2.49E-4 0.004547 1000 5001
Prediction using Eyjafjallajökull as a time-dependent covariate...
model {
for (i in 2:19) {
t[i] <- D[i] - D[i-1]
inactive[i] <- 1
erupt[i] <- 1
}
# Likelihood contributions from inactive periods
for (i in 2:19) {
inactive[i] ~ dbern(p.inactive[i])
p.inactive[i] <- exp(-mu*pow(t[i], r))
}
# Likelihood contributions at eruption times
for (i in 2:18) {
erupt[i] ~ dbern(q[i])
q[i] <- mu*r*pow(t[i], r-1)*exp(beta*eyja[i])
}
beta ~ dnorm(0, 0.7)
rel.risk <- exp(beta)
r ~ dunif(0, 10)
sigma ~ dgamma(0.001, 0.001)
mu <- 1/pow(sigma, r)
p.erupt.1 <- 1 - exp(mu*rel.risk*
(pow(92, r) - pow(92+1, r)))
p.erupt.50 <- 1 - exp(mu*rel.risk*
(pow(92, r) - pow(92+1, r)))*
exp(mu*
(pow(92+1, r) - pow(92+50, r)))
}
Data:
list(
D = c(1177, 1262 ,1311, 1357 , 1416 , 1440, 1450, 1500, 1550, 1580 , 1612, 1625, 1660, 1721, 1755, 1823, 1860, 1918, 2010),
eyja = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1))
Inits:
list(r=1, sigma=10, beta=0)
node mean sd MC error 2.5% median 97.5% start sample
beta 2.862 1.19 0.005721 0.5375 2.857 5.199 4001 216000
p.erupt.1 0.6546 0.2942 0.00139 0.09711 0.7047 1.0 4001 216000
p.erupt.50 0.983 0.0386 1.077E-4 0.8677 0.9981 1.0 4001 216000
r 2.112 0.4135 9.768E-4 1.365 2.09 2.983 4001 216000
rel.risk 35.6 61.91 0.2786 1.712 17.41 181.1 4001 216000
sigma 54.1 6.911 0.01656 41.69 53.68 68.95 4001 216000