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

model{

##prior, observer skill effect on detection, fixed
beta.p~dnorm(0,0.01)

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

### DP for occupancy
alpha~dgamma(a, b)  

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

##hyperparameters for detection component 
mu.a0 ~dnorm(0, 0.01)
sig.a0<-sqrt(1/tau.a0)
tau.a0~dgamma(0.01, 0.01)


##draw cluster level values from (mutivariate) DP base distributions
for (k in 1:K){
a0[k] ~ dnorm(mu.a0,tau.a0) #detection, univariate
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 stuff, cluster ID
	g.psi[i]~dcat(pi[1:K])

	a0[i] ~ dnorm(mu.a0,tau.a0) #detection, normal random effect


	for(j in 1:J){
	
	###detection stuff
	lp[i,j]<-a0[i]  + beta.p * OBS[j]  ##account for observer skill	
	p[i,j]<-pnorm(lp[i,j],0,1)

	lpsi[i,j]<-inprod(beta, VAR[j,]) + inprod(delta[g.psi[i],], VAR[j,])
	psi[i,j]<-pnorm(lpsi[i,j],0,1)
	z[i,j]~dbern(psi[i,j])
	p.eff[i,j]<-p[i,j]*z[i,j]
	y[i,j]~dbin(p.eff[i,j],K1[j])

	}

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

}


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

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

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

clusters.psi<-sum(cla[1:K])
}


















