INAR, adINAR, dpINAR & pyINAR

#library(devtools)
#install_github("heltongraziadei/HedINAR", auth_token = #"465593e8ea25032bfdda68c6fd58b49f5d32bd03")
library(HedINAR)

data("Pittsburgh")

year <- Pittsburgh$Year+(Pittsburgh$Month-1)/12
year[121]=2000
ts1 <- Pittsburgh$Area_12
ts2 <- Pittsburgh$Area_11

par(mfrow=c(1,1))
plot(year,ts1,type="b",ylim=range(ts1,ts2),ylab="Count")
lines(year,ts2,type="b",col=2)
legend("topright",legend=c("Area 12","Area 11"),col=1:2,lty=1)

boxplot(ts1[1:36],ts1[37:72],ts1[73:108],ts1[109:144],
        ts2[1:36],ts2[37:72],ts2[73:108],ts2[109:144],
        names=c("1-36","37-72","73-108","109-144",
                "1-36","37-72","73-108","109-144"),
        ylab="Count")
abline(v=4.5,lty=2)
text(2,33,"Area 12",col=2,cex=2)
text(6,33,"Area 11",col=2,cex=2)

inar1 <- inar(ts1)
inar2 <- inar(ts2)

names(inar1)

## [1] "time_series" "prior"       "burn_in"     "chain"       "est"        
## [6] "elapsed"

inar1$elapsed

##    user  system elapsed 
##   3.579   0.108   3.900

inar1$prior

## $a_alpha
## [1] 1
## 
## $b_alpha
## [1] 1
## 
## $a_lambda
## [1] 1
## 
## $b_lambda
## [1] 0.1

names(inar1$burn_in)

## [1] "length" "alpha"  "lambda"

names(inar1$chain)

## [1] "length" "alpha"  "lambda"

par(mfrow=c(2,3))
ts.plot(inar1$chain$alpha,xlab="Iteration",ylab="",main=expression(alpha))
acf(inar1$chain$alpha,main="")
hist(inar1$chain$alpha,xlab="",main=expression(alpha))
ts.plot(inar1$chain$lambda,xlab="Iteration",ylab="",main=expression(lambda))
acf(inar1$chain$lambda,main="")
hist(inar1$chain$lambda,xlab="",main=expression(lambda))

adinar1 <- adinar(ts1)
adinar2 <- adinar(ts2)

par(mfrow=c(2,2))
hist(adinar1$chain$alpha,xlab="",main=expression(alpha),prob=TRUE,xlim=c(0,0.6),ylim=c(0,10))
lines(density(inar1$chain$alpha),col=2)
hist(adinar1$chain$lambda,xlab="",main=expression(lambda),prob=TRUE,ylim=c(0,0.8),xlim=c(5,13))
lines(density(inar1$chain$lambda),col=2)
hist(adinar1$chain$w,xlab="",main=expression(w),prob=TRUE)
hist(1/adinar1$chain$theta,xlab="",main=expression(1/theta),
     breaks=seq(1,19,by=0.5),prob=TRUE,xlim=c(3,20))

par(mfrow=c(2,2))
hist(adinar2$chain$alpha,xlab="",main=expression(alpha),prob=TRUE,xlim=c(0,0.6),ylim=c(0,7))
lines(density(inar2$chain$alpha),col=2)
hist(adinar2$chain$lambda,xlab="",main=expression(lambda),prob=TRUE,ylim=c(0,2),xlim=c(1,4))
lines(density(inar2$chain$lambda),col=2)
hist(adinar2$chain$w,xlab="",main=expression(w),prob=TRUE)
hist(1/adinar2$chain$theta,xlab="",main=expression(1/theta),
     breaks=seq(1,19,by=0.25),prob=TRUE,xlim=c(0,10))

dpinar1 <- dpinar(ts1)
dpinar2 <- dpinar(ts2)

names(dpinar1)

## [1] "time_series" "prior"       "burn_in"     "chain"       "est"        
## [6] "elapsed"

names(dpinar1$chain)

## [1] "length"       "alpha"        "lambda"       "tau"         
## [5] "num_clusters"

par(mfrow=c(2,2))
hist(dpinar1$chain$alpha,prob=TRUE,xlab="",main=expression(alpha),xlim=c(0,0.5),ylim=c(0,10))
lines(density(inar1$chain$alpha),col=2)
lines(density(adinar1$chain$alpha),col=3)

hist(dpinar1$chain$tau,xlab="",main=expression(tau))

hist(dpinar1$chain$num_clusters,xlab="",main="clusters")

ts.plot(ts1[2:144],xlab="Time",ylab="")
lines(apply(dpinar1$chain$lambda,2,median),col=2)
legend("topright",legend=c("Time series","lambda (posterior median)"),lty=1,col=1:2)

par(mfrow=c(2,2))
hist(dpinar2$chain$alpha,prob=TRUE,xlab="",main=expression(alpha),xlim=c(0,0.5),ylim=c(0,10))
lines(density(inar2$chain$alpha),col=2)
lines(density(adinar2$chain$alpha),col=3)

hist(dpinar2$chain$tau,xlab="",main=expression(tau))

hist(dpinar2$chain$num_clusters,xlab="",main="clusters")

ts.plot(ts2[2:144],xlab="Time",ylab="")
lines(apply(dpinar2$chain$lambda,2,median),col=2)
legend("topright",legend=c("Time series","lambda (posterior median)"),lty=1,col=1:2)