nb10<-scan("nb10.txt")


qqnorm(nb10); qqline(nb10, col = 2)



# 1) MH=============================================

mu<-400
 n<-100
 u<-1/n*sum((nb10-mu)^2)
 s0<-60
 n0<-3
 f<-function(theta,n=n0,s=s0) 
     {
      (n/2)^(n/2)/gamma(n/2) * s^(n/2) * theta^(-n/2-1)*exp(-(n*s)/(2*theta))
     }

 likelihood<-function(sigma)
    {
      sigma^(-100/2) * exp(-sum((nb10-mu)^2)/(2*sigma))
    }

theta<-seq(from=0,to=200,length=1000)
plot(theta,f(theta,n-2,n*u/(n-2)), type='l',lty=1)
lines(theta,f(theta,n0,s0),lty=2)
lines(theta,f(theta,n+n0,(n0*s0+n*u)/(n0+n)),lty=3)

n0<-0.000001
plot(theta,f(theta,n-2,n*u/(n-2)), type='l',lty=1)
lines(theta,f(theta,n0,s0),lty=2)
lines(theta,f(theta,n+n0,(n0*s0+n*u)/(n0+n)),lty=3)

 rvscaleinvchisquare<-function(n,s)
 {
 x<-rchisq(1,n)
 x<-n*s/x
 return(x)
 }


 metropolis_hastings<-function(iter,nstar)
 {
  accept<-0
  met<-numeric(iter)
  last<-30
  for(i in 1:iter) {
  cand<-rvscaleinvchisquare(nstar,(nstar-2)/nstar * last)
  alpha<-(f(cand,(n0+100),(n0*s0+sum((nb10-mu)^2))/(n0+100))/f(last,(n0+100),(n0*s0+sum((nb10-mu)^2))/(n0+100)))*(f(last,nstar,((nstar-2)/nstar*cand))/f(cand,nstar,((nstar-2)/nstar*last)))
  if(runif(1)<min(alpha,1)){ 
  last<-cand 
  accept<-accept+1}
  met[i]<-last
 }
 list(sigma=met,accept=accept/iter)
 }

plot( metropolis_hastings(1000,20)$sigma, type="l")
plot( metropolis_hastings(1000,2.5)$sigma, type="l")
plot( metropolis_hastings(1000,200)$sigma, type="l")


# 2) MH (log) =============================================   

mu<-400
 s0<-0.01
 n0<-4.01

# for the proposal with normalizing constant

logf1<-function(theta,n=n0,s=s0) 
     {
     (n/2)*log(n/2) - lgamma(n/2) + (n/2)*log(s) + (-n/2-1)*log(theta) -(n*s)/(2*theta)
     }

# for the posterior without normalizing constant
logf<-function(theta,n=n0,s=s0) 
     {
     ( -n/2-1)*log(theta) -(n*s)/(2*theta)
     }

loglikelihood<-function(sigma)
    {
      (-100/2)*log(sigma) -sum((nb10-mu)^2)/(2*sigma) 
    }

 rvscaleinvchisquare<-function(n,s)
 {
 x<-rchisq(1,n)
 x<-n*s/x
 return(x)
 }


 metropolis_hastings2<-function(iter,nstar)
 {
  accept<-0
  met<-numeric(iter)
  last<-30
  for(i in 1:iter) {
  cand<-rvscaleinvchisquare(nstar,(nstar-2)/nstar * last)
  logalpha<-logf(cand,(n0+100),(n0*s0+sum((nb10-mu)^2))/(n0+100)) - logf(last,(n0+100),(n0*s0+sum((nb10-mu)^2))/(n0+100)) +logf1(last,nstar,((nstar-2)/nstar*cand)) -logf1(cand,nstar,((nstar-2)/nstar*last))
  if(  log(runif(1)) < logalpha ){ 
  last<-cand 
  accept<-accept+1}
  met[i]<-last
 }
 list(sigma=met,accept=accept/iter)
 }


 test3<-metropolis_hastings2( 1000, 500)
 test1<-metropolis_hastings2( 1000, 2.5)
 test2<-metropolis_hastings2( 1000, 20)

 par(mfrow=c(3,1))
 plot( test1$sigma, type='l')
 plot( test2$sigma, type='l')
 plot( test3$sigma, type='l')

# install InvChisquare {geoR}

metropolis_hastings3<-function(iter,nstar)
 {
  accept<-0
  met<-numeric(iter)
  last<-30
  for(i in 1:iter) {
  cand<-rinvchisq(1, nstar,(nstar-2)/nstar * last)
  logalpha<-dinvchisq(cand,(n0+100),(n0*s0+sum((nb10-mu)^2))/(n0+100), log=T) - dinvchisq(last,(n0+100),(n0*s0+sum((nb10-mu)^2))/(n0+100), log=T) +dinvchisq(last,nstar,((nstar-2)/nstar*cand), log=T) -dinvchisq(cand,nstar,((nstar-2)/nstar*last), log=T)
  if(  log(runif(1)) < logalpha ){ 
  last<-cand 
  accept<-accept+1}
  met[i]<-last
 }
 list(sigma=met,accept=accept/iter)
 }


# 3) MH (log) with normalizing constant ============================================= 

 mu<-400
 s0<-0.01
 n0<-4.01


logf<-function(theta,n=n0,s=s0) 
     {
     (n/2)*log(n/2) - lgamma(n/2) + (n/2)*log(s) + (-n/2-1)*log(theta) -(n*s)/(2*theta)
     }

loglikelihood<-function(sigma)
    {
      (-100/2)*log(sigma) -sum((nb10-mu)^2)/(2*sigma) 
    }

 rvscaleinvchisquare<-function(n,s)
 {
 x<-rchisq(1,n)
 x<-n*s/x
 return(x)
 }


 metropolis_hastings3<-function(iter,nstar)
 {
  accept<-0
  met<-numeric(iter)
  last<-30
  for(i in 1:iter) {
  cand<-rvscaleinvchisquare(nstar,(nstar-2)/nstar * last)
  logalpha<-logf(cand,(n0+100),(n0*s0+sum((nb10-mu)^2))/(n0+100)) - logf(last,(n0+100),(n0*s0+sum((nb10-mu)^2))/(n0+100)) +logf(last,nstar,((nstar-2)/nstar*cand)) -logf(cand,nstar,((nstar-2)/nstar*last))
  if(  log(runif(1)) < logalpha ){ 
  last<-cand 
  accept<-accept+1}
  met[i]<-last
 }
 list(sigma=met,accept=accept/iter)
 }


 test3<-metropolis_hastings3( 1000, 500)
 test1<-metropolis_hastings3( 1000, 2.5)
 test2<-metropolis_hastings3( 1000, 20)

 par(mfrow=c(3,1))
 plot( test1$sigma, type='l')
 plot( test2$sigma, type='l')
 plot( test3$sigma, type='l')


# 4) ergodic mean plot==============================================


erg.mean<-function(x){ # compute ergodic mean 
        n<-length(x)
        result<-cumsum(x)/cumsum(rep(1,n))
  }

em<-erg.mean(test2$sigma)

plot(em, type='l')


# 5) plots ========================================================


test2<-metropolis_hastings3( 1000, 20)

par( mfrow=c(3,1), xaxs='r', yaxs='r', bty='l' , cex=0.8)
iter<-1000
burnin<-100
index<-1:iter
index2<-(burnin+1):iter


plot(index, test2$sigma[index], type='l', ylab='Values of sigma', xlab='Iterations', main='(a) Trace Plot of sigma')
abline(v=100, col=2)



ergtheta1<-erg.mean( test2$sigma[index] )
ergtheta12<-erg.mean( test2$sigma[index2] )
ylims1<-range( c(ergtheta1,ergtheta12) )



step<-10
index3<-seq(1,iter,step)
index4<-seq(burnin+1,iter,step)


plot(index3, ergtheta1[index3], type='l', ylab='Values of sigma', xlab='Iterations', main='(b) Ergodic Mean Plot of sigma (step 10)', ylim=ylims1)
lines(index4, ergtheta12[index4-burnin], col=2, lty=2)
legend(400,55, lty=c(1,2), col=1:2,legend=c("Without removing Burnin Period", "Removing Burnin Period"))



acf(test2$sigma[index2], main='(c) Autocorrelations ')