#source("c:/Work07/ShortCourse07/Programs/Mastitis.txt",print.eval=TRUE)#
xdata <- c(0,0,1,1,2,2,2,2,2,2,4,4,4,5,5,5,5,5,5,6,6,8,8,8,9,9,9,10,10,12,12,13,13,13,13,18,18,19,19,19,19,20,20,22,22,22,23,25)
nx<-length(xdata)
nsim<-1000;
lambda<-array(2,dim=c(nsim,nx));beta<-array(5,dim=c(nsim,nx));
alpha<-.1;a<-1;b<-1;
for(i in 2:nsim){
for(j in 1:nx){
beta[i,j]<-1/rgamma(1,shape=alpha+a,scale=1/(lambda[i-1,j]+(1/b)));
lambda[i,j]<-rgamma(1,shape=xdata[j]+alpha,scale=1/(1+(1/beta[i,j])))
}
}
#**********************Plot Pre Processing***************************************
m<-250				#number of points to plot
step<-round(nsim/m) 		#nsim/m should be an integer!!!
lambda5<-lambda[,5][1];lambda15<-lambda[,15][1];beta15<-beta[,15][1];xplot<-1
for(i in 1:m) {lambda5<-c(lambda5,lambda[,5][i*step]);xplot<-c(xplot,i*step)}
for(i in 1:m) {lambda15<-c(lambda15,lambda[,15][i*step]);beta15<-c(beta15,beta[,15][i*step])}
den<-c(1:(m+1))
#**********************************************************************************
par(mfrow=c(2,3),mar=c(4.5,4.5,1.5,1.5))
hist(lambda5,main="Herd 5",freq=F,col="grey",breaks=25,xlim=c(0,8),xlab=expression(lambda[5]))
par(new=T)
plot(density(lambda5),main="",xlab="",ylab="",xaxt="n",yaxt="n",xlim=c(0,8),lwd=2)
hist(lambda15,main="Herd 15",freq=F,col="grey",breaks=30,xlim=c(0,12),xlab=expression(lambda[15]))
par(new=T)
plot(density(lambda15),main="",xlab="",ylab="",xaxt="n",yaxt="n",xlim=c(0,12),lwd=2)
hist(beta15,main="Herd 15",xlim=c(0,70),breaks=300,freq=F,col="grey",xlab=expression(beta[15]))
par(new=T)
plot(density(beta15),main="",xlab="",ylab="",xaxt="n",yaxt="n",xlim=c(0,70),lwd=2)
mean5<-cumsum(lambda5)/den
plot(xplot,mean5,type="l",col="black",xlim=c(0,nsim),ylab="Running Mean",xlab="Iteration",lwd=2)
mean15<-cumsum(lambda15)/den
plot(xplot,mean15,type="l",col="black",xlim=c(0,nsim),ylab="Running Mean",xlab="Iteration",lwd=2)
meanB15<-cumsum(beta15)/den
plot(xplot,meanB15,type="l",col="black",xlim=c(0,nsim),ylab="Running Mean",xlab="Iteration",lwd=2)