### We estimate effects of
### Rhodo, gap, year, transect etc based on logistic transformation.
### ...find.prob only for new seedlings.
### There are also functions for generating simulated data.
### Designed for the gap data, ie.40 quads by 6 years by 12 trans
### Random walk for true seedlings can move by 4.
### Data format:  c(36,12,6) for (quads,trans,years)
### Uses a separate lambda for the gap/canopy portion of each transect, i.e.,
### lambda is c(2,12,6-1) for c(c(can,gap),trans,years-1) since last year
### of data is not used.
### Updated 12/9/99.
#########################################################
############################################ FUNCTIONS
	
### FUNCTION 1. Simulates data:  first generate a data set, 
### then modify it to reflect probability of finding seedlings.

Simulate.gap<-function(n.quad,n.trans,n.years,surv.new,surv.old,f.prob){ #1
	# Simulates data of same form as gap plot data, ie
	# 12 transects of 32 quads x 6 yrs
	# f.prob prob of finding a seedling, same for all age classes
        # No treatment effects in simulated data
	
	#cat(paste(round(surv.new,2),round(surv.old,2),"\n", sep="\t"))
print(paste("Quadrats=", n.quad,"Surv.new=",surv.new,"Surv.old=",
	surv.old,"Find.prob=",f.prob,sep="  "))
	
new<-array(NA,c(n.quad,n.trans,n.years)) #  True New seedlings last year
old<-array(NA,c(n.quad,n.trans,n.years)) #  True Old seedlings last year
new.obs<-array(NA,c(n.quad,n.trans,n.years)) #  Observed New seedlings last year

for(j in 1:n.years){#3
for(i in 1:n.trans){#2
	
	rate<-4 #<-runif(1,min=2,max=6) #Each transect has different rate parameter
	new[,i,j]<-rpois(40,rate)
	
for(k in 1:n.quad){ #4
	
  if(j==1){old[k,i,j]<-rpois(1,2)}
  else{
  ### Generate true seedlings
  if(new[k,i,j-1]==0 && old[k,i,j-1]>0){old[k,i,j]<-rbinom(1,old[k,i,j-1],surv.old)}
  if(new[k,i,j-1]>0 && old[k,i,j-1]==0){old[k,i,j]<-rbinom(1,new[k,i,j-1],surv.new)}
  if(new[k,i,j-1]>0 && old[k,i,j-1]>0){old[k,i,j]<-rbinom(1,new[k,i,j-1],surv.new)+
     rbinom(1,old[k,i,j-1],surv.old)}
  if(new[k,i,j-1]==0 && old[k,i,j-1]==0){old[k,i,j]<-0}
  } #end else
  
  ### Generate observed seedlings for new seedlings only
  if(new[k,i,j]>0) {new.obs[k,i,j]<-rbinom(1,new[k,i,j],f.prob)}
  	else {new.obs[k,i,j]<-0}
  }#4
  }#2
  }#3
return(new,old,new.obs)
} #1 end function

#################### generate data with treatment effects
Simulate.gap.beta<-function(f.prob,X.mat,n.quad=32,n.trans=12,n.years=6){ #1
	# Simulates data of same form as gap plot data, ie
	# 12 transects of 32 quads x 6 yrs
	# f.prob prob of finding a seedling, same for all age classes
	# Treatment effects are included.
	#cat(paste(round(surv.new,2),round(surv.old,2),"\n", sep="\t"))
print(paste("Quadrats=", n.quad,"Surv.new=",surv.new,"Surv.old=",
	surv.old,"Find.prob=",f.prob,sep="  "))

##### lambda's
lam.n<-array(NA,c(2,n.trans,n.years))
lam.n[1,,]<-runif(n.trans*n.years,min=2,max=6) # lambda's for canopy
lam.n[2,,2:n.years]<-lam.n[1,,2:n.years]+rnorm(n.trans*(n.years-1), mean=1, sd=.2) # lambda's for gap
	lam.n[2,,1]<-lam.n[1,,1]+rnorm(n.trans, mean=0, sd=.2) # lambda's for gap in year 1

#### Beta's
beta<-matrix(NA,nrow=1,ncol=22)
	dimnames(beta)<-list(1,c("Bo","Bg","Br","Bgr","Bt1",
	"Bt2","Bt3","Bt4","Bt5","Bt6","Bt7","Bt8","Bt9","Bt10","Bt11","Bt12",
	"By1","By2","By3","By4","By5","Bd"))
beta[1]<--2.4
beta[2]<-0.5
beta[3]<--1
beta[4]<-0
beta[22]<-3
beta[5:16]<-rnorm(12,0,0.4)
beta[17:21]<-rnorm(5,0,0.4)

ETAstar.new<-ETA.new<-array(NA,c(32,12,5))
for(i in 1:12){
	for(j in 1:5){ # We can only use n.yr-1 or 5 years of data
		ETA.new[,i,j]<-X.mat[,i,j,]%*%beta[-22]
		}}		# B[22] is Bd
ETA.old<-ETA.new+beta[22]

Sn<-exp(ETA.new)/(1+exp(ETA.new))
Sd<-exp(ETA.old)/(1+exp(ETA.old))
	
new<-array(NA,c(n.quad,n.trans,n.years)) #  True New seedlings last year
old<-array(NA,c(n.quad,n.trans,n.years)) #  True Old seedlings last year
new.obs<-array(NA,c(n.quad,n.trans,n.years)) #  Observed New seedlings last year

for(i in 1:n.trans){#2
	for(j in 1:n.years){#3

	new[c(1:8,25:32),i,j]<-rpois(16,lam.n[1,i,j])
	new[9:24,i,j]<-rpois(16,lam.n[2,i,j])
	
for(k in 1:n.quad){ #4
	
  if(j==1){old[k,i,j]<-rpois(1,2)}
  else{
  ### Generate true seedlings
  if(new[k,i,j-1]==0 && old[k,i,j-1]>0){old[k,i,j]<-rbinom(1,old[k,i,j-1],Sd[k,i,j-1])}
  if(new[k,i,j-1]>0 && old[k,i,j-1]==0){old[k,i,j]<-rbinom(1,new[k,i,j-1],Sn[k,i,j-1])}
  if(new[k,i,j-1]>0 && old[k,i,j-1]>0){old[k,i,j]<-rbinom(1,new[k,i,j-1],Sn[k,i,j-1])+
     rbinom(1,old[k,i,j-1],Sd[k,i,j-1])}
  if(new[k,i,j-1]==0 && old[k,i,j-1]==0){old[k,i,j]<-0}
  } #end else
  
  ### Generate observed seedlings for new seedlings only
  if(new[k,i,j]>0) {new.obs[k,i,j]<-rbinom(1,new[k,i,j],f.prob)}
  	else {new.obs[k,i,j]<-0}
  }#4
  }#2
  }#3
return(new,old,new.obs,beta)
} #1 end function

########################################
### FUNCTION to create design matrix

X.mat<-function(){
# I don't include Bd in design array...it's dealt with separately

	X<-array(0,c(32,12,5,21)) # Doesn't include year 6 or Bd
	dimnames(X)<-list(1:32,1:12,1:5,c("Bo","Bg","Br","Bgr","Bt1",
	"Bt2","Bt3","Bt4","Bt5","Bt6","Bt7","Bt8","Bt9","Bt10","Bt11","Bt12",
	"By1","By2","By3","By4","By5"))
X[,,,1]<-1
X[9:24,,2:4,2]<-1 # 1~gap:  there was a gap in years 2-4 but not 1
X[,c(4:6,10:12),,3]<-1; X[,c(1:3,7:9),,3]<-0 # 1~ Rhodo
X[,,,4]<-X[,,,3]*X[,,,2]
for(i in 1:12) { X[,i,,i+4]<-1}
for(i in 1:5){ X[,,i,i+16]<-1}

return(X)
}

########################################
### FUNCTION called by Surv.lik
	
Surv.lik.c<-function(data.array, llik.array, llik.sum, n.qd, n.tr,n.yr, Sn, Sd){
	.C("likelihood_beta",as.integer(data.array),
	llik.array,llik.sum,n.qd,n.tr,n.yr,as.numeric(Sn),as.numeric(Sd),
	classes=c("integer","numeric","numeric","integer","integer",
	"integer","numeric","numeric"))[[2]]}

########################################
### FUNCTION 2 Calculates log likelihood of all
### quadrats given true seedling numbers and survival probs.
		
Surv.lik_function(new,old,surv.new,surv.old){
	#Surv.lik(N,D,Snstar,Sd)  called from Beta.samp, N.samp
	# recall that N ~dim(n.qd,n.tr,n.yr-1) but D~dim(n.qd,n,tr,n.yr)

	n.quads<-dim(new)[1]; n.trans<-dim(new)[2]; n.years<-dim(new)[3]
	data<-array(NA,c(n.quads,n.years,3,n.trans))
for(i in 1:n.trans){
  data[,,1,i]<-new[,i,];  data[,,2,i]<-old[,i,1:n.years]
  data[,,3,i]<-old[,i,2:(n.years+1)]
	}
	# Data structure allows each vector of length 3 ie data[1,1,]
	# to have new.yr0,old.yr0, and old.yr1
 
trash2<-array(9999.0,c(32,12,5))
trash3<-numeric(1)

	# print("calling Surv.lik.c") 
llik.quads<-Surv.lik.c(data,trash2,trash3,n.quads,n.trans,n.years,surv.new,surv.old)
	# print("returning from Surv.lik.c")
		
	llik.quads[llik.quads==9999]<-NA
llik.sum<-sum(llik.quads,na.rm=T)

return(llik.sum,llik.quads)
 } # End Function


#################################################################
############################################### Sampling functions

F.samp<-function(n,N,alpha.prior,beta.prior){
	# function call does not pass last year of data in n
	# which is consistent with dim of imputed N
	
	obs.sum<-numeric(1); true.sum<-numeric(1)
obs.sum<-sum(n); true.sum<-sum(N)
alpha<-(obs.sum+alpha.prior)
beta<-(true.sum-obs.sum+beta.prior)
out<-rbeta(1,alpha,beta)
	
#cat(paste("f.sim",round(out,2),"true.sum",round(true.sum,2),
#	"Alpha",round(alpha,2),"Beta",round(beta,2),"\n",sep="\t"))
	
return(out)
}

##########
Lam.samp<-function(n.qd,n.tr,n.yr,N,alpha.prior,beta.prior){
	# Lam.samp is partly hardwired.	Doesn't use final year of data.
lam<-array(NA,c(2,n.tr,(n.yr-1)))	
for(i in 1:n.tr){
	for(j in 1:(n.yr-1)){
		# rgamma(1,alpha)/beta
	lam[1,i,j]<-rgamma(1,sum(N[c(1:8,25:32),i,j])+alpha.prior)/(16+beta.prior)
		# Canopy
	lam[2,i,j]<-rgamma(1,sum(N[9:24,i,j])+alpha.prior)/(16+beta.prior)
		# Gap
	}}

return(lam)
}

##########
Rt.samp<-function(n,df,mu,std){
	# samples from unidimensional t
m<-length(mu) #ie length=1
y<-rt(n,df)
y<-y*(std)
y<-mu+y
return(y)
}

#########
Dt.samp<-function(x,df=100,mu,std){
	# evaluates density of univariate t
m<-length(mu)
tmp<-(x-mu)/std
y<-dt(tmp,df)/std
return(y)}

######### FUNCTION called by Beta.samp

Design.c<-function(eta.array, bstar, b, X, n.qd, n.tr,n.yr, n.beta){
	.C("designmat_fcn",eta.array,bstar,b,as.integer(X),
	n.qd, n.tr, n.yr, n.beta,
	classes=c("numeric","numeric","numeric","integer","integer",
	"integer","integer","integer"))[[1]]}

##########
Beta.samp<-function(N,D,B,k,X,mu.prior,sigma2.prior,std=0.2,scale=0.5){
	#Sn<-Bo+Br+Bg+Bgr+Bt[1:12]+By[1:6]
	#Sd<-Bo+Br+Bg+Bgr+Bt[1:12]+By[1:6]+Bd

beta.star<-Rt.samp(1,100,B[k],std*scale)
Bstar<-B; Bstar[k]<-beta.star

ETA<-array(9999.0,c(32,12,5,2)) # I use 9999 to indicate an NA

tmp<-Design.c(ETA,Bstar[-22],B[-22],X,32,12,5,21)

ETA.new<-tmp[,,,1];  ETAstar.new<-tmp[,,,2]
ETAstar.old<-ETAstar.new+Bstar[22];  ETA.old<-ETA.new+B[22]
			# B[22] is Bd

Snstar<-exp(ETAstar.new)/(1+exp(ETAstar.new))
Sdstar<-exp(ETAstar.old)/(1+exp(ETAstar.old))

Sn<-exp(ETA.new)/(1+exp(ETA.new))
Sd<-exp(ETA.old)/(1+exp(ETA.old))


part1a<-Surv.lik(N,D,Snstar,Sdstar)$llik.sum
part1b<-log(dnorm(Bstar[k], mu.prior, sd=sqrt(sigma2.prior)))

part2a<-Surv.lik(N,D,Sn,Sd)$llik.sum
part2b<-log(dnorm(B[k], mu.prior, sd=sqrt(sigma2.prior)))

Pi.dif<-exp(part1a+part1b-part2a-part2b)	# because I was using log likelihoods
alpha<-min(1,Pi.dif)

if(alpha>=runif(1)){
	out<-list(B=Bstar[k],Sn=Snstar,Sd=Sdstar)} 
else{
	out<-list(B=B[k],Sn=Sn,Sd=Sd)} 
	#eq to (runif(1)<=alpha)
	
	#cat(paste("Snstar",round(Snstar,2),"Pi.dif",round(Pi.dif,2),
	#	"Alpha",round(alpha,2),"Sn",round(out,2),"\n",sep="\t"))

return(out)
}

##########
Sigma2.samp<-function(N,D,B,mu,alpha.prior,beta.prior){
	# B should only be the relevant portion of the B vector
	alpha<-beta<-out<-numeric()
n<-length(B) # equals n-1 because last year already not considered
alpha<-alpha.prior+0.5*n
beta<-beta.prior + 0.5*t(B-mu)%*%(B-mu)

out<-rgamma(1,alpha)/beta
out<-1/out			#InvGamma(alpha,beta)
return(out)
}

##########
Q.fun<-function(z,Z,zstar,n.quads,n.trans,n.years){
# Q.fun calculates the edge correction for random walk (max step 2)
# in N.samp,D.samp
# Z~true seedlings (imputed data),
# z~observed seedlings (yr 0 does not matter here)
# zstar~candidate point

	Q<-array(0,c(n.quads,n.trans,n.years))
	
	Q[Z>z+3]<-0 #Shouldn't need this if Q default is 0
	Q[Z==zstar]<-0 # Shouldn't need this if Q default is 0
	Q[(Z==z & zstar==(z+1))]<-log(3/4)
	Q[Z==z & zstar==(z+2)]<-log(3/5)
	Q[Z==(z+1) & zstar==z]<-log(4/3)
	Q[Z==(z+1) & zstar==(z+2)]<-log(4/5)
	Q[Z==(z+1) & zstar==(z+3)]<-log(4/5)
	Q[Z==(z+2) & zstar==z]<-log(5/3)
	Q[Z==(z+2) & zstar==(z+1)]<-log(5/4)
	Q[Z==(z+3) & zstar==(z+1)]<-log(5/4)

return(Q)
}

########################################
### FUNCTION called by N.samp.  It calls a c function.

Q.func.c<-function(z,Z,Zstar,Q,n.quads,n.trans,n.years){
	.C("Qfunction",as.integer(z),as.integer(Z),as.integer(Zstar),
	Q,n.quads,n.trans,n.years,classes=c("integer","integer","integer",
	"numeric","integer","integer","integer"))[[4]]}


##########
N.samp<-function(n,N,D,lam.n,f.prob,Sn,Sd){
	# N.samp(n,N,D,lam.n[,,,i],fn[i],tmp$Sn,tmp$Sd)
	# recall that N ~dim(n.qd,n.yr-1) but n,D ~dim(n.qd,n.yr)

	n.quads<-dim(N)[1]; n.trans<-dim(N)[2]  ; n.years<-dim(N)[3]
	n.quad.total<-n.quads*n.trans*n.years
	
dif<-N-n[,,1:n.years]

nstar<-N+ifelse(dif>=4,sample(-4:4,n.quad.total,replace=T),
	ifelse(dif==3,sample(-3:4,n.quad.total,replace=T),
	ifelse(dif==2,sample(-2:4,n.quad.total,replace=T),
	ifelse(dif==1,sample(-1:4,n.quad.total,replace=T),	
	sample(0:4,n.quad.total,replace=T)))))

part1a<-Surv.lik(nstar,D,Sn,Sd)$llik.quads  # array(40,12,5)	
part1b<-array(log(dbinom(n[,,1:n.years],nstar,f.prob)),
	c(n.quads,n.trans,n.years))

part2a<-Surv.lik(N,D,Sn,Sd)$llik.quads	
part2b<-array(log(dbinom(n[,,1:n.years],N,f.prob)),
	c(n.quads,n.trans,n.years))

part1c<-part2c<-array(NA,c(n.quads,n.trans,n.years))

	for(i in 1:n.trans){
		for(j in 1:n.years){
part1c[c(1:8,25:32),i,j]<-log(dpois(nstar[c(1:8,25:32),i,j],lam.n[1,i,j])) # Canopy
part1c[9:24,i,j]<-log(dpois(nstar[9:24,i,j],lam.n[2,i,j])) # Gap
part2c[c(1:8,25:32),i,j]<-log(dpois(N[c(1:8,25:32),i,j],lam.n[1,i,j])) # Canopy
part2c[9:24,i,j]<-log(dpois(N[9:24,i,j],lam.n[2,i,j])) # Gap
			}}
	# OLD: part2c<-array(log(dpois(N,lam.n)),c(n.quads,n.trans,n.years))
	
Q<-array(0.0,c(32,12,5))
part3Q<-Q.func.c(n[,,1:n.years],N,nstar,Q,n.quads,n.trans,n.years)	#New c function
	
part1a[is.na(part1a)]<-0	;part2a[is.na(part2a)]<-0	
part1b[is.na(part1b)]<-0	;part2b[is.na(part2b)]<-0
		 
part4.sum<-part1a+part1b+part1c-part2a-part2b-part2c+part3Q
	
	tmp2<-array(NA,c(n.quads,n.trans,n.years,2)); tmp2[,,,1]<-1 
	tmp2[,,,2]<-exp(part4.sum)
alpha<-pmin(tmp2[,,,1],tmp2[,,,2])
	N.old<-N
N<-ifelse(runif(n.quad.total)<=alpha,nstar,N) # Check this!
	
#bug.check<-min((N[,1]+D[,1])-D[,2],(N[,2]+D[,2])-D[,3],(N[,3]+D[,3])-D[,4])

return(N)
}

################################################### Boa.fcn
### Combines chains into a single matrix
	
Boa.fcn<-function(fn,N.total,Sigma2,B)
	{
n<-length(fn)

mcmc.boa<-cbind(fn,N.total,Sigma2,B)
	
dimnames(mcmc.boa)<-list(1:n,c("fn","N.total","S2.t","S2.y","Bo","Bg","Br","Bgr","Bt1",
	"Bt2","Bt3","Bt4","Bt5","Bt6","Bt7","Bt8","Bt9","Bt10","Bt11","Bt12",
	"By1","By2","By3","By4","By5","Bd"))

return(mcmc.boa)
	}

################################################ Main.fcn
### Main function creates the iterations
Main.fcn<-function(mcmc.init,lam.init,X,N,n,D,n.sample){
	priter<-100; counter<-1 # keep track of progess
	
####### Creating objects ###
fn<-numeric(n.sample+1); lam.n<-array(NA,c(2,n.tr,(n.yr-1),n.sample+1))

B<-matrix(NA,nrow=(n.sample+1),ncol=22)
	dimnames(B)<-list(1:(n.sample+1),c("Bo","Bg","Br","Bgr","Bt1",
	"Bt2","Bt3","Bt4","Bt5","Bt6","Bt7","Bt8","Bt9","Bt10","Bt11","Bt12",
	"By1","By2","By3","By4","By5","Bd"))
	
Sigma2<-matrix(NA,nrow=(n.sample+1),ncol=2)
	dimnames(Sigma2)<-list(1:(n.sample+1),c("Bt","By"))

N.total<-numeric(n.sample+1) # for debugging
		
####### Initialization

mu.prior.years<-rep(0,5); mu.prior.trans<-rep(0,12)
	lam.n[,,,1]<-lam.init
	B[1,]<-mcmc.init[5:26]
	fn[1]<-mcmc.init[1]; N.total[1]<-mcmc.init[2]
	Sigma2[1,]<-mcmc.init[3:4]
	# N remains as N from previous run

####### ITERATION LOOP	

for(i in 2:(n.sample+1)){ #1 

	##### SAMPLING Fn: Gibbs step
fn[i]<-F.samp(n[,,1:(n.yr-1)],N,1,1) # only uses n.yr-1 years of data
	#print("finished F.samp")

	##### SAMPLING lam.n:  Gibbs steps
lam.n[,,,i]<-Lam.samp(n.qd,n.tr,n.yr,N,0.06,0.01)
	#Lam.samp<-function(n.qd,n.tr,n.yr,N,alpha.prior,beta.prior)
	#print("finished Lam.samp")	

	##### SAMPLING Beta's using random walk
	for(j in 1:22){
		if(j<=4 || j==22){
			tmp<-Beta.samp(N,D,B[i-1,],j,X,0,1000)} # Non Random effects
		else{
			if(j<=16 && j>=5){
				tmp<-Beta.samp(N,D,B[i-1,],j,X,0,Sigma2[i-1,1])}
					# transects random effects
			else{
				tmp<-Beta.samp(N,D,B[i-1,],j,X,0,Sigma2[i-1,2])}
					# year random effects
			}	
B[i,j]<-tmp$B
	#Beta.samp<-(N,D,B,k,X,mu.prior,sigma2.prior,std=0.2,scale=0.5)
	#returns out<-list(B=Bstar,Sn=Snstar,Sd=Sdstar)} 
	#cat(paste("finished Beta.samp",j),"\n",sep="\t")
	}

	##### SAMPLING Sigma's
	##### dimnames(Sigma2)<-list(1:(n.sample+1),c("Bt","By"))
Sigma2[i,1]<-Sigma2.samp(N,D,B[i,5:16],mu.prior.trans,0.001,0.001) 
	# sigma for TRANSECTS:  Sigma2.samp(N,D,B,mu,alpha.prior,beta.prior)  
	#print("finished Sigma2.samp ")

Sigma2[i,2]<-Sigma2.samp(N,D,B[i,17:21],mu.prior.years,0.001,0.001)  
	# sigma for YEARS:  Sigma2.samp(N,D,B,mu,alpha.prior,beta.prior)
	#print("finished Sigma2.samp ")
	
	##### Sampling N using random walk
N<-N.samp(n,N,D,lam.n[,,,i],fn[i],tmp$Sn, tmp$Sd)
	# N.samp(n,N,D,B,lam.n,f.prob,Sn,Sd)
	#print("finished N.samp")

	N.total[i]<-sum(N) #part of debugging code

if(priter*counter==(i+1)){print(paste("sample #",i+1,sep=" "));counter<-counter+1}


} #1 end of iter loop

	mcmc<-Boa.fcn(fn,N.total,Sigma2,B)
	mcmc<-mcmc[-1,]; lam.n<-lam.n[,,,-1]
return(mcmc,lam.n,N)
} # End of Main.fcn

#################################################################
################################################ BODY OF PROGRAM
#################################################################

Gap.beta2.fcn<-function(n,D,n.samples,z2,mcmc.out){
### z2==1~initialize with N.old, mcmc.old,lam.old
### z2==0~generate initial values

### Creating the data
	z<-0
	if(z==1){
	print("Generating data...")
### SET DATA GENERATION PARAMTERS!! ###
n.quad<-32; n.trans<-12; n.years<-6; surv.new<-0.1; surv.old<-0.7; f.prob<-0.8
	X<-X.mat() 
data<-Simulate.gap.beta(f.prob,X,n.quad=32,n.trans=12,n.years=6)
	#data<-Simulate.gap(n.quad,n.trans,n.years,surv.new,surv.old,f.prob)

n<-data$new.obs
Nt<-data$new  ### Can't refer to this as N because N is the imputed data
D<-data$old
}

################################################ Initial values
	print("Creating initial values and storage objects...")

### Setting Universal values
	#n.samples<-10               ### SET NUMBER OF ITERATIONS!!! ##
		iter.main<-500 # number of iterations per main call
	n.main.calls<-n.samples/iter.main
	n.qd<-dim(n)[1] # row dimension
	n.tr<-dim(n)[2] # colum dimension
	n.yr<-dim(n)[3] # 3rd dimension # dim(N) is n.yr-1
	
####### Creating objects ###
lam.n<-array(NA,c(2,n.tr,(n.yr-1),n.samples+1))

mcmc<-matrix(NA,nrow=(n.samples+1),ncol=26)
dimnames(mcmc)<-list(1:(n.samples+1),c("fn","N.total","S2.t","S2.y","Bo","Bg","Br","Bgr","Bt1",
	"Bt2","Bt3","Bt4","Bt5","Bt6","Bt7","Bt8","Bt9","Bt10","Bt11","Bt12",
	"By1","By2","By3","By4","By5","Bd"))
################################### INITIAL VALUES ###
#z2<-0  if z2==0 then initialize with new parameter values
if(z2==0){
	tmp<-array(NA,c(n.qd,n.tr,(n.yr-1))) # Holds initial values for imputed 
		# number of new new seedlings	

	for(i in 1:(n.yr-1)){tmp[,,i]<-D[,,i+1]-(n[,,i]+D[,,i])}
	tmp<-ifelse(tmp<0,0,tmp) ; N<-n[,,1:(n.yr-1)]+tmp 
		# Note:	N has dim(n.qd,n.tr,n.yr-1)
	#N<-Nt[,,1:(n.yr-1)] # can only be used with simulated data

	for(i in 1:n.tr){
		for(j in 1:(n.yr-1)){
	lam.n[1,i,j,1]<-mean(n[1:16,i,j])
	lam.n[2,i,j,1]<-mean(n[1:16,i,j])
	}}

	mcmc[1,5:26]<-0;	mcmc[1,1]<-.75
	mcmc[1,3:4]<-1000; mcmc[1,2]<-sum(N)
	print("Intialized with new values")	
	}
else{
	mcmc.old<-mcmc.out$mcmc; lam.old<-mcmc.out$lam.n
	tmp1<-dim(mcmc.old)[1]
	lam.n[,,,1]<-lam.old[,,,tmp1]
	mcmc[1,5:26]<-mcmc.old[tmp1,5:26] #B
	mcmc[1,1]<-mcmc.old[tmp1,1] #fn
	mcmc[1,2]<-mcmc.old[tmp1,2] #N.total
	mcmc[1,3:4]<-mcmc.old[tmp1,3:4] #Sigma2
	N<-mcmc.out$N # N remains as N from previous run
	print("Intialized with old values")
}
				
################################### ITERATIONS ###
print(paste(n.samples,"samples to be generated",sep=" "))

for(r in 1:n.main.calls){
	d<-((r-1)*iter.main+1)
		cat(paste("Beginning iteration",d),sep="  ","\n")
	out<-Main.fcn(mcmc[d,],lam.n[,,,d],X,N,n,D,iter.main)
		#Main.fnc(mcmc.init,lam.init,X,N,n.sample){
	mcmc[(d+1):(d+iter.main),]<-out$mcmc
	lam.n[,,,c(d+1):(d+iter.main)]<-out$lam.n
	N<-out$N
}
	cat("Iterations completed, exiting successfully...","\n",sep="\t")
	
return(mcmc,lam.n,N)
} # End Gap.beta2.fcn