####################################################
#Importance sampling: Gamma

#version I - know normalizing constants Zp, Zq
alpha<-7.8
beta<-2

n<-100000
a<-floor(alpha)
b<-beta*a/alpha 
U<-matrix(rexp(a*n),a,n)
Y<-apply(U,2,'sum')/b

exp_Y1<-mean(Y*Y^(alpha-a)*exp(-(beta-b)*Y)*(gamma(a)/gamma(alpha))*(beta^alpha/b^a))

exp_Y1
alpha/beta

#version II dont know normalizations
num_Y<-mean(Y*Y^(alpha-a)*exp(-(beta-b)*Y))
den_Y<-mean(Y^(alpha-a)*exp(-(beta-b)*Y))
exp_Y2<-num_Y/den_Y

exp_Y2
alpha/beta

den_Y
(gamma(alpha)/gamma(a))*(b^a/beta^alpha)

#################################################################
#see the estimator break down as alpha moves away from a

alpha<-seq(from=1,to=15,length.out=m)
beta<-2

n<-1000
a<-7
b<-1
U<-matrix(rexp(a*n),a,n)
Y<-apply(U,2,'sum')/b
m<-20
ex<-rep(NA,m)
#calulate the mean of f(X)=1 (so E(f)=1) in Gamma(alpha,beta) using Gamma(a,b)
for (i in 1:m) {
	ex[i]<-mean(Y^(alpha[i]-a)*exp(-(beta-b)*Y)
			*(gamma(a)/gamma(alpha[i]))*(beta^alpha[i]/b^a))
}
plot(alpha,ex)
lines(alpha,alpha^0)

####################################################
#IS - normal example, method I, Pr(Z>3) Z~N(0,1)
#simple approach (of lectures) without discussion of tilting (as in notes)
#note I chose x=4 so we could check using the naive method
#for x=6 say the advantage is even more dramatic

sigma<-1
mu<-0
x0<-4

n<-1000000
Z<-rnorm(n)
(plain_est<-mean(Z>x0))
(plain_sd<-sd(Z>x0)/sqrt(n))
 
Y<-x0+sigma*rnorm(n)
wpq<-(Y>x0)*exp(-(Y-mu)^2/(2*sigma^2)+(Y-x0)^2/(2*sigma^2))
(IS_est<-mean( wpq ))
(IS_std<-sd( wpq )/sqrt(n))

####################################################
#IS - normal example, method I, Pr(Z>3) Z~N(0,1)
#uses tilting - this is the version in lecture notes

sigma<-1
mu<-0
x<-4

n<-1000000
Z<-rnorm(n)
(plain_est<-mean(Z>x))
(plain_sd<-sd(Z>x)/sqrt(n))

t<-(x-mu)/sigma^2
Y<-mu+t*sigma^2+sigma*rnorm(n)
wpq<-(Y>x)*exp(-(Y-mu)^2/(2*sigma^2)+(Y-mu-t*sigma^2)^2/(2*sigma^2))
(IS_est<-mean( wpq ))
(IS_std<-sd( wpq )/sqrt(n))

####################################################
#MCMC simulate X_t according to a Poisson dbn of mean lambda=3.
lambda<-3

n<-10000
X<-rep(NA,n)
X[1]<-0
for (t in 1:(n-1)) {
	x<-X[t]
	y<-x+2*(runif(1)<=0.5)-1
	if (y>x) {
		a<-lambda/(x+1)
	}
	if (y<x & y>=0) {
		a<-x/lambda
	}
	if (y<0) {
		a<-0
	}
	U<-runif(1)
	if (U<=a) {
		X[t+1]<-y
	}
	else {
		X[t+1]<-x
	}
}

mean(X)
var(X)
windows(10,5); par(mfrow=c(1,2))
plot(X[1:200],type="l")
hist(X,-0.5:max(X+0.5),freq=F); lines(0:max(X),dpois(0:max(X),lambda))

############################################################
#rejection for X~p, p=[1:m]/sum(1:m).
m<-30
p<-(1:m)/sum(1:m)
q<-rep(1/m,m)
M<-max(p/q)

n<-1000
U<-runif(n)
Y<-ceiling(m*runif(n))
ap<-(p[Y]/(M*q[Y]))
X<-Y[U<ap]

par(mfcol=c(2,2))
plot(X,type='l',ann=F)
hist(X,-1:m+0.5,freq=F,ann=F); lines(1:m,pt/sum(pt))

#
#MCMC simulate X_t according to p=[1:m]/sum(1:m).
m<-30
pt<-1:m

n<-10000
X<-rep(NA,n)
X[1]<-1
for (t in 1:(n-1)) {
	x<-X[t]
	y<-ceiling(m*runif(1))
	a<-min(1,y/x)
	U<-runif(1)
	if (U<=a) {
		X[t+1]<-y
	}
	else {
		X[t+1]<-x
	}
}

plot(X[1:200],type="l",ann=F)
hist(X,-1:m+0.5,freq=F,ann=F); lines(1:m,pt/sum(pt),ann=F)
####################################################

####################################################
#MCMC simulate Normal mu=5,sigma^2=1
mu<-5
Ssq<-1

n<-10000
X<-rep(NA,n)
X[1]<-0
a<-1

for (t in 1:(n-1)) {
	x<-X[t]
	y<-x+(2*runif(1)-1)*a
	U<-runif(1)
	MHR<-exp( ((x-mu)^2-(y-mu)^2)/(2*Ssq) )
	alpha<-min(1,MHR)
	if (U<=alpha) {
		X[t+1]<-y
	}
	else {
		X[t+1]<-x
	}
}

par(mfcol=c(2,3))
plot(X[1:1000],type='l')
h<-hist(X,freq=F); 
lines(h$mids,dnorm(h$mids,mu,sqrt(Ssq)))
mean(X)
var(X)

####################################################
#MCMC simulate bivariate Normal
mu<-c(1,1)
(S<-matrix(c(3,-2,-2,3),2,2))
(Si<-solve(S))

n<-10000
X<-matrix(c(NA,NA),2,n)
X[,1]<-c(0,0)
a<-3

for (t in 1:(n-1)) {
	x<-X[,t]
	y<-x+(2*runif(2)-1)*a
	U<-runif(1)
	MHR<-exp( -t(y-mu)%*%Si%*%(y-mu)/2+t(x-mu)%*%Si%*%(x-mu)/2 )
	alpha<-min(1,MHR)
	if (U<=alpha) {
		X[,t+1]<-y
	}
	else {
		X[,t+1]<-x
	}
}

plot(t(X[,1:10]),type='l',xlim=c(-8,8),ylim=c(-8,8))
plot(t(X[,1:20]),type='l',xlim=c(-8,8),ylim=c(-8,8))
plot(t(X[,1:100]),type='l',xlim=c(-8,8),ylim=c(-8,8))
plot(t(X[,1:10000]),xlim=c(-8,8),ylim=c(-8,8))


apply(X,1,mean)
cov(t(X))

#see notes
plot(t(X[,1:100]),type='l',xlab='X_1(t)',ylab='X_2(t)',cex.lab=1.5,xlim=c(-6,6),ylim=c(-4,6))
plot(t(X[,1:2000]),xlab='X_1(t)',ylab='X_2(t)',cex.lab=1.5,xlim=c(-6,6),ylim=c(-4,6),col='red')


####################################################
#MCMC simulate bivariate Normal - conditioned |X-(3,3)|<1
mu<-c(1,1)
(S<-matrix(c(3,-2,-2,3),2,2))
(Si<-solve(S))

n<-10000
X<-matrix(c(NA,NA),2,n)
X[,1]<-c(3,3) #X_0 at center of allowed region
a<-1

for (t in 1:(n-1)) {
	x<-X[,t]
	y<-x+(2*runif(2)-1)*a
	U<-runif(1)
	MHR<-exp( -t(y-mu)%*%Si%*%(y-mu)/2+t(x-mu)%*%Si%*%(x-mu)/2 )
	alpha<-min(1,MHR)
	if (U<=alpha & sum((y-c(3,3))^2)<1 ) {
		X[,t+1]<-y
	}
	else {
		X[,t+1]<-x
	}
}
 
points(t(X[,1:2000]),pch=20,xlab='X_1(t)',ylab='X_2(t)',cex.lab=1.5,xlim=c(-6,6),ylim=c(-4,6))
 
theta<-seq(from=0,to=2*pi,length.out=30)
v<-3+cos(theta)
w<-3+sin(theta)
lines(v,w,col='blue',lwd=1)

##############################################################
#Ising model
n<-64
theta<-0.8

#X0 is random binary vector
X<-matrix(as.numeric(runif(n^2)>0.5),n,n)
dx<-sum(abs(diff(X))+abs(diff(t(X))))
N<-1000000
for (j in 1:N) {
    i<-1+floor(runif(1)*n^2) #choose a component to update
    Y<-X
    Y[i]<-1-X[i]
    dy<-sum(abs(diff(Y))+abs(diff(t(Y)))) #d(y)
    if (runif(1)<exp(theta*(dx-dy))) {
        X<-Y
        dx<-dy       #d(y) is saved as d(x) if y is accepted
    }
    if (floor(j/1000)==j/1000) {image(X,col=gray(0:255/255),axes=F)}
}