# CHAPTER 8: R COMMANDS
# 08/01/2007

# IMAGE SEGMENTATION PART

###################################################################################

source("#8.R")

gisi1=isinggibbs(1000,100,0.3)
gisi2=isinggibbs(1000,100,0.4)
gisi3=isinggibbs(1000,100,0.5)
gisi4=isinggibbs(1000,100,0.6)
gisi5=isinggibbs(1000,100,0.7)
gisi6=isinggibbs(1000,100,0.8)
gisi7=isinggibbs(1000,100,0.9)
gisi8=isinggibbs(1000,100,1)
gisi9=isinggibbs(1000,100,1.1)
gisi10=isinggibbs(1000,100,1.2)
par(mfrow=c(5,2),mar=c(2,2,2,2))
image(1:100,1:100,t(gisi1),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(gisi6),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(gisi2),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(gisi7),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(gisi3),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(gisi8),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(gisi4),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(gisi9),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(gisi5),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(gisi10),col=gray(c(1,0)),xlab="",ylab="")
dev.copy2eps(file="gising.eps")

###################################################################################

hmisi1=isinghm(1000,100,0.3)
hmisi2=isinghm(1000,100,0.4)
hmisi3=isinghm(1000,100,0.5)
hmisi4=isinghm(1000,100,0.6)
hmisi5=isinghm(1000,100,0.7)
hmisi6=isinghm(1000,100,0.8)
hmisi7=isinghm(1000,100,0.9)
hmisi8=isinghm(1000,100,1)
hmisi9=isinghm(1000,100,1.1)
hmisi10=isinghm(1000,100,1.2)
par(mfrow=c(5,2),mar=c(2,2,2,2))
image(1:100,1:100,t(hmisi1),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(hmisi6),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(hmisi2),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(hmisi7),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(hmisi3),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(hmisi8),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(hmisi4),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(hmisi9),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(hmisi5),col=gray(c(1,0)),xlab="",ylab="")
image(1:100,1:100,t(hmisi10),col=gray(c(1,0)),xlab="",ylab="")
dev.copy2eps(file="hmising.eps")

###################################################################################

gpotts1=pottsgibbs(1000,100,0.3)
gpotts2=pottsgibbs(1000,100,0.4)
gpotts3=pottsgibbs(1000,100,0.5)
gpotts4=pottsgibbs(1000,100,0.6)
gpotts5=pottsgibbs(1000,100,0.7)
gpotts6=pottsgibbs(1000,100,0.8)
gpotts7=pottsgibbs(1000,100,0.9)
gpotts8=pottsgibbs(1000,100,1.0)
gpotts9=pottsgibbs(1000,100,1.1)
gpotts10=pottsgibbs(1000,100,1.2)
par(mfrow=c(5,2),mar=c(2,2,2,2))
image(1:100,1:100,t(gpotts1),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(gpotts6),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(gpotts2),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(gpotts7),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(gpotts3),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(gpotts8),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(gpotts4),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(gpotts9),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(gpotts5),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(gpotts10),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
dev.copy2eps(file="gpotts.eps")

###################################################################################

hmpotts1=pottshm(1000,100,0.3)
hmpotts2=pottshm(1000,100,0.4)
hmpotts3=pottshm(1000,100,0.5)
hmpotts4=pottshm(1000,100,0.6)
hmpotts5=pottshm(1000,100,0.7)
hmpotts6=pottshm(1000,100,0.8)
hmpotts7=pottshm(1000,100,0.9)
hmpotts8=pottshm(1000,100,1.0)
hmpotts9=pottshm(1000,100,1.1)
hmpotts10=pottshm(1000,100,1.2)
par(mfrow=c(5,2),mar=c(2,2,2,2))
image(1:100,1:100,t(hmpotts1),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(hmpotts6),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(hmpotts2),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(hmpotts7),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(hmpotts3),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(hmpotts8),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(hmpotts4),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(hmpotts9),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(hmpotts5),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
image(1:100,1:100,t(hmpotts10),col=gray(c(1,0.75,0.25,0)),xlab="",ylab="")
dev.copy2eps(file="hmpotts.eps")

###################################################################################

matu=matrix(0,2000,21)
for (i in 1:21) 
{
matu[,i]=pottshm6(2000,100,0.1*i-0.1)
print(i)
}

matu=matrix(0,2000,21)
for (i in 1:21)
{
matu[,i]=pottshm6(2000,100,0.1*i-0.1)
print(i)
}

m=rep(0,21)
for (i in 1:21) m[i]=mean(matu[1001:2000,i])
lrcst=approxfun(seq(0,2,0.1),m)

plot(seq(0,2,0.1),m,xlab="",ylab="")
curve(lrcst,0,2,add=T)
dev.copy2eps(file="pathsampling.eps")

integrate(lrcst,0.3,0.4)$value

###################################################################################

lm3=read.table("Menteith")
image(1:100,1:100,lm3,col=gray(256:1/256),xlab="",ylab="")
dev.copy2eps(file="menteith.eps")

###################################################################################

truncnorm=function(n,mu,tau2,a,b)
{
qnorm(pnorm(b,mu,sqrt(tau2))-runif(n)*(pnorm(b,mu,sqrt(tau2))-pnorm(a,mu,sqrt(tau2))),mu,sqrt(tau2))
}

###################################################################################

titus=reconstruct1(2000,lm3,m)

affect=function(u)
{
order(u)[6]
}

aff=apply(titus$xcum,1,affect)
aff=t(matrix(aff,100,100))
par(mfrow=c(2,1))
image(1:100,1:100,aff,col=gray(6:1/6),xlab="",ylab="")
image(1:100,1:100,lm3,col=gray(256:1/256),xlab="",ylab="")
dev.copy2eps(file="recons.eps")
dev.off()

par(mfrow=c(3,2),mar=c(2,2,2,2))
plot(titus$mu[101:2000,1],type="l")
plot(titus$mu[101:2000,2],type="l")
plot(titus$mu[101:2000,3],type="l")
plot(titus$mu[101:2000,4],type="l")
plot(titus$mu[101:2000,5],type="l")
plot(titus$mu[101:2000,6],type="l")
dev.copy2eps(file="mut.eps")
dev.off()

par(mfrow=c(3,2),mar=c(2,2,2,2))
hist(titus$mu[101:2000,1],nclass=100,main="")
hist(titus$mu[101:2000,2],nclass=100,main="")
hist(titus$mu[101:2000,3],nclass=100,main="")
hist(titus$mu[101:2000,4],nclass=100,main="")
hist(titus$mu[101:2000,5],nclass=100,main="")
hist(titus$mu[101:2000,6],nclass=100,main="")
dev.copy2eps(file="muh.eps")
dev.off()

par(mfrow=c(2,2),mar=c(2,2,2,2))
plot(titus$sigma2[101:2000],type="l")
hist(titus$sigma2[101:2000],nclass=50,main="")
plot(titus$beta[101:2000],type="l")
hist(titus$beta[101:2000],nclass=50,main="")
dev.copy2eps(file="betesig.eps")
dev.off()

###################################################################################