parameter estimation of rossler equation by using rjags

rjags_rossler.R

library(rjags)

alpha<-1/10.
beta<-1/10.
c<-6

N <- 500
dt<-0.1
sig<-1e-4
itenum<-500000

x<-c()
y<-c()
z<-c()

x[1]<-5
y[1]<-1
z[1]<-0.
  
for (t in 2:N) {
  x[t] <- x[t-1]+dt*(-y[t-1]-z[t-1])+rnorm(1,0,1e-3)
  y[t] <- y[t-1]+dt*(x[t-1]+alpha*y[t-1])+rnorm(1,0,1e-3)
  z[t] <- z[t-1]+dt*(beta+z[t-1]*(x[t-1]-c))+rnorm(1,0,1e-3)
}

#plot(x,type='l')

set.seed(21)
model <- "
var                                                                                                                                                                                    
N, x[N],y[N],z[N], alpha,beta,c,dt,sig,out
model {                     
x[1] ~ dnorm(10.0, 1.0E-4)   
y[1] ~ dnorm(10.0, 1.0E-4)   
z[1] ~ dnorm(0.0, 1.0E-4)   
for (t in 2:N) {            
  x[t] ~ dnorm(x[t-1]+dt*(-y[t-1]-z[t-1]),sig)
  y[t] ~ dnorm(y[t-1]+dt*(x[t-1]+alpha*y[t-1]),sig)
  z[t] ~ dnorm(z[t-1]+dt*(beta+z[t-1]*(x[t-1]-c)),sig)
}
c ~ dnorm(0.0, 1.0E-3)
out <- x[N]
}"

write(model, "tmp.bug")


init1 <- list(c=5)
init2 <- list(c=7)
init3 <- list(c=10)
inits<-list(init1, init2, init2)

m <- jags.model("tmp.bug",
                data = list(N=N,z=z,alpha=alpha,beta=beta,dt=dt,sig=sig),
                inits = inits,
                n.chains = length(inits), n.adapt = itenum/10)

post <- coda.samples(m,
                     variable.names = c("c", "x","out"),
                     n.iter = itenum, thin = 20)

post.sum <- summary(post)
post.var <- varnames(post)

isplot<-T

#ggmcmc
library(ggmcmc)
ggr<-ggs(post)
ggs_density(ggr,family="out")
ggs_traceplot(ggr,family="out")


## plot                                                                                                                                                                              
if(isplot){
  x.ind <- grep("x\\[", post.var)
  x.mean <- post.sum$statistics[x.ind, "Mean"]
  x.025 <- post.sum$quantiles[x.ind, "2.5%"]
  x.975 <- post.sum$quantiles[x.ind, "97.5%"]  
  
plot(1:N, x, type = "b", las = 1,xlab = "time", ylab = "n")
lines(1:N, x.mean, col = 2)
lines(1:N, x.025, col = 2, lty = 2)
lines(1:N, x.975, col = 2, lty = 2)
}