pdf(file="chap4-gibbs.pdf", paper="a4")
##################################################
#Suppose we want to simulate from a centered bivariate normal (μ=0) with correlation ρ and variances both 1 and we know that X|y∼N(ρy,1−ρ2) 
##################################################
### example taken from http://web.stanford.edu/class/stats366/Gibbs.html ### 

######Function######
gibbsNimate <- function(NSimul=50, rho=0.5, x0=1, y0=1)
{
  ##NSimul=Number of simulations
  ##rho is the correlation
  ##x0 and y0 are the starting values
  ## plot the 95% contour of the theoretical normal
  ##
  plot(ellipse(rho), type="l",
   main="Bivariate Normal", col="green")
      dev.hold()
  ## initialization
  x <- rep(0, NSimul)
  y <- rep(0, NSimul)
  x[1]=x0
  y[1]=y0
  ## Gibbs sampler and plot points
#
  for(b in 2:NSimul) {
     plot(ellipse(rho), type="l",
               main="Bivariate Normal")
  points(x[1:(b-1)], y[1:(b-1)])
  dev.hold()
    ## sample the x direction conditional on y
  x[b] <- rnorm(1, rho*y[b-1], sqrt(1-rho^2))
  lines(c(x[b-1], x[b]), c(y[b-1], y[b-1]))
  ## sample the y direction conditional on x
  y[b] <- rnorm(1, rho*x[b], sqrt(1-rho^2))
  lines(c(x[b], x[b]), c(y[b-1], y[b]))
  ani.pause()
  points(x[b], y[b],col="red")
  }
  ## return the draws
  return(cbind(x,y))
}


#######
require(ellipse)
require(animation)
oopt = ani.options(interval = 2)
gibbsNimate(20,rho=0.5)


################## nuclear pump failure data ################################
failures=c( 5,1,5,14,3,19,1,1,4,22)
times=c(94.32, 15.72, 62.88, 125.76, 5.24, 31.44, 1.05, 1.05, 2.1, 10.48)
pump=data.frame(failures,times)
print(pump)

nx=length(failures)
Nsim=5e4

lambda=array(NA, c(Nsim, nx))
beta=numeric(Nsim)

lambda[1,]=rep(1,nx)
beta[1]=1
alpha=1.8
gamma=.01
delta=1
for(i in 2:Nsim){
        lambda[i,]=lambda[i-1,]
        beta[i]=beta[i-1]
    for(j in 1:nx){
        lambda[i,j] = rgamma(1, sh=failures[j]+alpha, ra=times[j]+beta[i])
    }
        beta[i] = rgamma(1, sh=nx*alpha + gamma, ra=delta + sum(lambda[i,]))
}

plot.ts(lambda[,])
plot.ts(beta)
op=par(mfrow=c(5,2), mar=c(2, 2, 1, 1))
for(j in 1:10)
hist(lambda[,j], breaks='FD')
par(op)

hist(beta, breaks='FD')

colMeans(lambda)
apply(lambda, 2, mean)

acf(lambda[,])
pacf(lambda[,1])



apply(lambda^2, 2, mean)

dev.off()

############################## Ising model on a 'line' ############################## 

sum.neighbours = function(S, i){
    K = length(S)
    if( (i==1) ){
        my.sum= S[i+1]
    } else if( i==K ){
        my.sum= S[i-1]
    } else {
        my.sum= S[i+1] + S[i-1]
    }
    return(my.sum)
}


simulate.ising.gibbs = function(S0, K, H, J, T){
    S=S0
#
for(t in 2:T){
    ## print   progress
    if((t %% floor(1+T/200))==0){ cat("\n",round(100*t/T,1), "%--", sep="") }
    for(i in 1:K){
            sum.neigh.tmp = sum.neighbours(S,i)
            if(is.na(sum.neigh.tmp)) stop()
            p =  1/ (1 +
                                                                 exp(2*(J*sum.neigh.tmp
                                                                        + H)))
            if(is.na(p)) stop()
            S[i] = rbinom(1,1, p)*2 - 1
    }
}
return(S)
}


K=20
J=-5
H=5
T=1e3
S0=rbinom(K, 1, .5)*2-1
#S=simulate.ising.gibbs(S0, K, H, J, T)

pdf(file="chap4-ising-1d.pdf", paper="a4")
image(as.matrix(S0))
#
op=par(mfrow=c(2,2))
system.time(
for(J in c(-20, -1, 1, 20)){
#
S=simulate.ising.gibbs(S0, K=K, J=J, H=H, T=T)
#
#image(S[T-100,,], main=paste("J = ", J, "H = ", H,", T=",  T-1))
image(as.matrix(S), main=paste("J = ", J, "H = ", H, ",T=",  T), breaks=c(-2,0,2), col=c(1,8))
#
}
)
par(op)
dev.off()

############################## Ising Model ##############################
sum.neighbours = function(S, i,j){
    K = dim(S)[2]
    if( (i==1) ){
        if( j==1 ){
        my.sum= S[i+1,j] + S[i, j+1]
        }else if (j == K) {
        my.sum= S[i+1,j] + S[i, j-1]
        } else { ## j not on the edge
        my.sum= S[i+1,j] + S[i, j-1] + S[i, j+1]
        }
    } else if( i==K ){
        if( j==1 ){
        my.sum= S[i-1,j] + S[i, j+1]
        }else if (j == K) {
        my.sum= S[i-1,j] + S[i, j-1]
        } else { ## j not on the edge
        my.sum= S[i-1,j] + S[i, j-1] + S[i, j+1]
        }
    } else {
        if( j==1 ){
        my.sum= S[i+1,j] + S[i-1,j] + S[i, j+1]
        }else if (j == K) {
        my.sum= S[i+1,j] + S[i-1,j] + S[i, j-1] 
        } else { ## j not on the edge
        my.sum= S[i+1,j] + S[i-1,j] + S[i, j-1] + S[i, j+1]
        }
    }
    return(my.sum)
}

simulate.ising.gibbs = function(S0, K, H, J, T){
    S=S0
#
for(t in 2:T){
    ## print   progress
    if((t %% floor(1+T/200))==0){ cat("\n",round(100*t/T,1), "%--", sep="") }
    for(i in 1:K){
        for(j in 1:K){
            sum.neigh.tmp = sum.neighbours(S,i,j)
            if(is.na(sum.neigh.tmp)) stop()
            p = 1 / (1 + exp(2*(J*sum.neigh.tmp + H)))
            if(is.na(p)) stop()
            S[i,j] = rbinom(1,1, p)*2 - 1
        }
    }
}
return(S)
}



old.simulate.ising.gibbs = function(K, H, J, T){
#
S= array(NA, c(T, K, K)) ## variable to store the markov chain
S[1,,]= matrix(-1 + 2*rbinom(K^2, 1, .5), K, K) ## initial value of the MC
for(t in 2:T){
    ## print   progress
    if((t %% floor(1+t/200))==0){ cat("\n",round(100*t/T,1), "%--", sep="") }
    S[t,,] = S[t-1,,]
    for(i in 1:K){
        for(j in 1:K){
            sum.neigh.tmp = sum.neighbours(S[t,,],i,j)
            if(is.na(sum.neigh.tmp)) stop()
            p = ( 1 ) / (1 + exp(2*(J*sum.neigh.tmp + H)))
            if(is.na(p)) stop()
            S[t,i,j] = rbinom(1,1, p)*2 - 1
        }
    }
}
return(S)
}





pdf(file="chap4-ising.pdf", paper="a4")

K=100
H=0
p0=.5
S0=matrix(rbinom(K^2, 1, p0)*2-1, K, K)
T=1000
image(S0,  breaks=c(-2,0,2), col=c(1,8))
#
op=par(mfrow=c(2,2))
time.tot1=system.time(
for(J in c(-10, -1, 0,1)){
#
#
S=simulate.ising.gibbs(S0, K=K, J=J, H=H, T=T)
#
image(S[,], main=paste("J = ", J, "H = ", H, ",T=",  T), breaks=c(-2,0,2), col=c(1,8) )
#
}
)
par(op)

K=100
H=0
p0=.5
S0=matrix(rbinom(K^2, 1, p0)*2-1, K, K)
image(S0,  breaks=c(-2,0,2), col=c(1,8))
#
op=par(mfrow=c(2,2))
time.tot2=system.time(
for(J in c(-1, -.5, -.3, -.1)){
#
S=simulate.ising.gibbs(S0, K=K, J=J, H=H, T=T)
#
image(S[,], main=paste("J = ", J, "H = ", H, ",T=",  T), breaks=c(-2,0,2), col=c(1,8) )
#
}
)
par(op)

dev.off()