#
# This program calculates marginal likelihoods for two different software reliability models and 
# computes the log Bayes factor.
#
# You should open the file containing the data and mark and copy this.  Then after the first command, 
# paste the data directly.
#
musa <- read.table(stdin())
count <- musa$V1
x <- musa$V2
#
# First the independent exponential model
#
a <- 1
b <- 1
apost <- a+length(x)
bpost <- b+sum(x)
fx <- exp(log(apost)+apost*log(bpost)-(apost+1)*log(bpost+xgrid))
likind <- lgamma(a+length(x))-lgamma(a)+a*log(b)-(a+length(x))*log(b+sum(x))
etexp <- bpost/(apost-1)
experror <- mean((x-etexp)^2)
#
# Now a Jelinski Moranda model.
#
lambda <- 200
ajm <- 1
bjm <- 1
iters <- 10000
n <- rep(0,iters)
theta <- rep(0,iters)
lengthx <- length(x)
n[1] <- lengthx
minn <- lengthx
sumx <- sum(x)
sumix <- sum(c(1:lengthx)*x)
#
# Initial Gibbs sampler run.
#
theta[1] <- rgamma(1,lengthx,(n[1]+1)*sumx-sumix)
et <- 1/(theta[1]*((n[1]+1)*sumx-sumix))
for (i in 2:iters){
  n[i] <- lengthx+rpois(1,exp(log(lambda)-theta[i-1]*sumx))
  theta[i] <- rgamma(1,ajm+lengthx,bjm+(n[i]+1)*sumx-sumix)
  et <- et+1/(theta[i]*((n[i]+1)*sumx-sumix))
}
#
# Plot of posterior densities.
#
hist(n,probability=TRUE,main='',xlab='N',ylab='p',breaks=c(length(x):max(n)))
et <- et/iters
jmerror <- mean((et-x)^2)
#
# Marginal likelihood evaluation.
#
nmode <- as.numeric(names(sort(-table(n))[1]))
etheta <- mean(theta)
logpri <- dpois(nmode,lambda,log=TRUE)+dgamma(etheta,a,b,log=TRUE)
loglik <- lengthx*log(etheta)+lgamma(nmode+1)-lgamma(nmode-lengthx+1)-etheta*((nmode+1)*sumx-sumix)
n[1] <- n[iters]
theta[1] <- theta[iters]
logftheta <- dgamma(etheta,ajm+lengthx,bjm+(n[1]+1)*sumx-sumix,log=TRUE)
#
# Second sampler run.
#
for (i in 2:iters){
  n[i] <- lengthx+rpois(1,exp(log(lambda)-theta[i-1]*sumx))
  theta[i] <- rgamma(1,ajm+lengthx,bjm+(n[i]+1)*sumx-sumix)
  logftheta <- logftheta+dgamma(etheta,ajm+lengthx,bjm+(n[i]+1)*sumx-sumix,log=TRUE)
}
logftheta <-logftheta/iters
logfn <- dpois(nmode-lengthx,exp(log(lambda)-etheta*sumx),log=TRUE)
likjm <- logpri+loglik-logftheta-logfn
#
# Now Bayes factor.
#
bayesfactor <- likjm-likind
experror
jmerror