data{
for (xx in 1:npar){
mu[xx]<-0  ##mean for G prior
mu.b[xx]<-0
}
}

model{

##prior for detection
p~dbeta(1,1)

##prior, g-prior parameter
omega~dt(0,1,1)T(0.001,) #scaled T-distr. (half), with phi=1 and df=1 (Johnson example script)

alpha~dgamma(a,b)T(0.001,)  ###Devin Johnson's prior

##stick-breaking stuff for parameter alpha
	pp[1] ~ dbeta(1,alpha)T(0.001,0.999)
	pi[1] <- pp[1]
	for (k in 2:(K.max-1)){
		pp[k] ~ dbeta(1,alpha)T(0.001,0.999)
		pi[k] <- pp[k]*(1-pp[k-1])*pi[k-1]/pp[k-1]
	}
	ps<-sum(pi[1:(K.max-1)])
	pi[K.max] <- 1-ps


###draw cluster level values from mutivariate DP base distributions
for (k in 1:K.max){
delta[k,1:npar] ~ dmnorm.vcov(mu[1:npar],Omega.mat[1:npar,1:npar])
}

for (xx in 1:npar){
for (yy in 1:npar){
Omega.mat[xx,yy]<-omega^2 * R[xx,yy] ##R = (H'H)^-1; data; g-prior 
}}


##fixed parameters
beta[1:npar]~dmnorm.vcov(mu.b[1:npar],R.beta[1:npar,1:npar])


for (i in 1:n){

	### occupancy, cluster ID
	g[i]~dcat(pi[1:K.max])

	for(j in 1:J){
	lpsi[i,j]<-inprod(delta[g[i],],X[j,]) + inprod(beta,X[j,]) #frst term: species level effects, second term: fixed (avg) effects  
	psi[i,j]<-pnorm(lpsi[i,j],0,1)  #probit link
	z[i,j]~dbern(psi[i,j])
	p.eff[i,j]<-p*z[i,j]
	y[i,j]~dbin(p.eff[i,j],K)

	}

Nocc[i]<-sum(z[i,1:J])

}


###table specs per cluster to get at number of cluster

for (k in 1:K.max){
	for(i in 1:n){
	SCb[k,i]<-equals(g[i], k)
}}

for (k in 1:K.max){
clb[k]<-step(sum(SCb[k,1:n])-1)  ###step() has to be negative to return 0
}

clusters<-sum(clb[1:K.max])

}


















