#source("c:/Work07/ShortCourse07/Programs/Chinook.txt",print.eval=TRUE)#
#Estimates parameters with Gibbs sampling
#Data for Year 2
source("c:/Work07/ShortCourse07/Programs/Subroutines.txt",print.eval=TRUE)#
data<-read.table("c:/Work07/ShortCourse07/Programs/ChinookData.txt",sep = "",header=T)
Length<-data[,1]
Year<-data[,2]
Age<-data[,3]
#--------------------Select Year and Order Observaions--------------------
#--------------------So missing data are together-------------------------
Length<-Length[Year==2]
Age<-Age[Year==2]
n<-length(Age);nm<-sum(is.na(Age))		#count missing data
oa<-order(Age)
Length<-Length[oa];Age<-Age[oa]		#missing data at end
#--------------------Prior Parameters-------------------------------------
m0<-c(400,600,800,1000,1100)			#Length prior mean
tau2=c(3,6,6,6,5)				#Length prior variance
a<-3;b<-1/100					#Inverted Gamma parameters, sets mean = 50
#--------------------Initial Values---------------------------------------
nsim<-10000		#number of Gibbs iterations
p<-array(0,dim=c(nsim,5))
nj<-(n/5)*c(1,1,1,1,1)
p[1,]<-rdirichlet(1,(n/5)*c(1,1,1,1,1))
mu<-m0
sigma2<-1/rgamma(5,a,scale=b)
#--------------------Gibbs Iterations-------------------------------------
for(i in 1:nsim)
{ print(i)
  p[i,]<-rdirichlet(1,nj+1)		#Generate p
#--------------------Fill in Missing Age----------------------------------
  for(j in (n-nm+1):n)
  {
    pAge<-p[i,]*dnorm(Length[j],mean=mu,sd=sqrt(sigma2))
    pAge<-pAge/sum(pAge)
    Age[j]<-sum(c(1:5)*rmultinom(1,1,pAge))	#multinomial on 1-5 scale
  }
#--------------------Summary Statistics-----------------------------------
  nj<-tapply(Length,Age,length)
  xbar<-tapply(Length,Age,mean)
  s2<-tapply(Length,Age,var)
#--------------------Generate Parameters-----------------------------------
  sigma2<-1/rgamma(5,(nj/2)+a,rate=nj*s2/2+1/b)
  w<-nj*tau2/(sigma2+nj*tau2)
  mu<-rnorm(5,mean=w*xbar+(1-w)*m0,sd=sqrt(tau2*(1-w)))
}
#--------------------End of Gibbs Iterations-------------------------------
#--------------------Plot Output-------------------------------
par(mfrow=c(2,3))
for(i in 1:5){hist(p[,i],col="grey",xlab="",main=i+2)}
pmean<-apply(p,2,cumsum)/c(1:nsim)
for(j in 1:5){plot(c(1:nsim),pmean[,j],type="l",ylim=c(0,.5));par(new=T)}