# Sample R program to simulate and analyze base model.
# The base model is 0-censored normal distribution.  Some values are also missing.

# Setup
library('coda')

rm(list=ls())
setwd("~/Data/ScytheMCMC/examples")

# Set numeric missing value that matches base_model.h missing definition.
MISSING <- "-9999"

# Number of observations and missing observations
N <- 1000
N.missing <- 0
censor <- FALSE;
# Sample from 0-censored normal with mean=0.25, sd = 1.
X <- rnorm(N, .25, 1)
if(censor) {
  X[X<0] = 0 # left censor at 0.
}
if (N.missing>0) {
  X[1:N.missing] <- NA # Add missing data  
}
hist(X)

# Save table in Scythe format (dense matrix with spaces)
write.table(X, 'X.txt', quote=FALSE, na=MISSING, row.names=FALSE, col.names=FALSE, sep=" ")
system("./base_model --config-file=base_model.ini --chains=1 --sample-size=5000 --burnin=5000 --thin=0 --random-seed=123456 --out-file=out.txt")

result <- read.csv('out.txt')
resultmc <- as.mcmc(result)
summary(resultmc)
plot(resultmc)
autocorr.diag(resultmc)
rejectionRate(resultmc)
effectiveSize(resultmc)

hist(result$mean)
abline(v=.25)