#A Gibbs sampler for the BLNE model without individual heterogeneity:


#Number of samples in chain
L=220000

#Length of burn-in
burn=20000

#Set random walk length for \omega_U
omega_U=5

#Original Data
y=c(3,1,2,1,1,3,2,1,3,1,3,5,1,1,4,5,3,6,2,2,5,5,5)
k=7
T=45

#Augmented data with M-nstar individuals with y=0
M=80
nstar=length(y)
Aug=M-nstar
y=c(y,rep(0,Aug))

#Define function "DataAugBin" for MCMC algorithm

DataAugBin=function(L,y,k,T,M,omega_U)	{

	mu=rep(0,L)
	q=rep(0,M)
	nstar=sum(y!=0)
	q[1:nstar]=1
	psi=rep(0,L)
	n=rep(0,L)
	N=rep(0,L)
	U=rep(0,L)
	mu[1]=runif(1)
	psi[1]=runif(1)
	n[1]=nstar
	maxU=1000 #set maximum value for U[1]
	U[1]=sample(seq(ceiling((T/k)),maxU),1)
	N[1]=U[1]+n[1]

	for (i in 2:L)	{
		psi[i]=rbeta(1,(1+sum(q)),(1+M-sum(q)))
		pq=(1-mu[i-1])^k*psi[i]/((1-mu[i-1])^k*psi[i]+(1-psi[i]))
		q[(nstar+1):M]=rbinom((M-nstar),1,pq)
		n[i]=sum(q)
		mu[i]=rbeta(1,1+sum(y)+T,1+k*(n[i]+U[i-1])-sum(y)-T)
		
		#M-H step on U using uniform random walk
		Ustar=U[i-1]+sample(seq(-omega_U,omega_U,1),1)
		if (Ustar*k<T)	{R=0} else {
			fUstar<-dbinom(T,k*(Ustar),mu[i])*1/Ustar
			fU<-dbinom(T,k*(U[i-1]),mu[i])*1/U[i-1]
			R<-fUstar/fU
		}
		U[i]<-ifelse(runif(1) < min(1,R),Ustar,U[i-1])
		N[i]=U[i]+n[i]


	}

	out=list(N,mu,psi,n)
	return(out)
}


#A hybrid Metropolis-within-Gibbs algorithm for the BLNE model with individual heterogeneity:

#Number of samples in chain
L=220000

#Length of burn-in
burn=20000

#Original Data
y=c(3,1,2,1,1,3,2,1,3,1,3,5,1,1,4,5,3,6,2,2,5,5,5)
k=7
U=45

#Augmented data with M-nstar individuals with y=0
M=80
nstar=length(y)
Aug=M-nstar
y=c(y,rep(0,Aug))

#Specify values for hyperprior parameters
tau=0.25
A=0.5
B=0.5

#Specify random walk step lengths
omega_U=5
omega_theta=2
omega_beta=0.2
omega_sigma=0.15

#Specify r for multiple-try MH step
r=4

#Define function "DataAugLogit" for MCMC algorithm

DataAugLogit=function(L,tau,A,B,y,k,T,M,omega_U,omega_theta,omega_beta,omega_sigma,r)    {

	theta=matrix(0,nrow=L,ncol=M)
	beta=rep(0,L)
	sigma2=rep(0,L)
	q=rep(0,M)
	nstar=sum(y!=0)
	Q=matrix(0,nrow=L,ncol=(M-nstar))
	q[1:nstar]=1
	Q[1,]=q[(nstar+1):M]
	psi=rep(0,L)
	mu=rep(0,L)
	n=rep(0,L)
	U=rep(0,L)
	N=rep(0,L)
	beta[1]=rnorm(1,0,sqrt(tau))
	sigma2[1]=rgamma(1,A,B)
	theta[1,]=rnorm(M,beta[1],(sqrt(sigma2[1])))
	psi[1]=rbeta(1,(1+sum(q)),(1+M-sum(q)))
	n[1]=sum(q)
	mu[1]=1/(1+exp(-beta[1]))
	maxU=1000 #set maximum value for U[1]
	U[1]=sample(seq(ceiling((T/k)),maxU),1)
	N[1]=U[1]+n[1]

	for (i in 2:L)  {

		#Draw from univariate conditional distributions of psi and q
		psi[i]=rbeta(1,(1+sum(q)),(1+M-sum(q)))
		pq=(1-1/(1+exp(-theta[i-1,])))^k*psi[i]/((1-1/(1+exp(-theta[i-1,])))^k*psi[i]+(1-psi[i]))
		q[(nstar+1):M]=rbinom((M-nstar),1,pq[(nstar+1):M])
		Q[i,]=q[(nstar+1):M]
		n[i]=sum(q)

		#MH step for theta_s using uniform random walk
		thetastar=theta[i-1,]+runif(M,-omega_theta,omega_theta)
		pstar=1/(1+exp(-thetastar))
		p=1/(1+exp(-theta[i-1,]))
		fthetastar=(pstar*q)^y*(1-pstar*q)^(k-y)*dnorm(thetastar,beta[i-1],(sqrt(sigma2[i-1])))
		ftheta=(p*q)^y*(1-p*q)^(k-y)*dnorm(theta[i-1,],beta[i-1],(sqrt(sigma2[i-1])))
		R=rep(0,M)
		R[fthetastar>=ftheta]=1
		R[fthetastar<ftheta]=(fthetastar/ftheta)[fthetastar<ftheta]
		probthetastar=runif(M)<R
		theta[i,][probthetastar==TRUE]=thetastar[probthetastar==TRUE]
		theta[i,][probthetastar==FALSE]=theta[i-1,][probthetastar==FALSE]

		#Draw from univariate conditional distribution of beta
		betastar=beta[i-1]+runif(1,-omega_beta,omega_beta)
		mustar=1/(1+exp(-betastar))
		fbetastar=prod(dnorm(theta[i,],betastar,(sqrt(sigma2[i-1]))))*dnorm(betastar,0,sqrt(tau))*dbinom(T,k*(U[i-1]),mustar)
		fbeta=prod(dnorm(theta[i,],beta[i-1],(sqrt(sigma2[i-1]))))*dnorm(beta[i-1],0,sqrt(tau))*dbinom(T,k*(U[i-1]),mu[i-1])
		R=fbetastar/fbeta
		if (runif(1) < min(1,R))	{
			beta[i]=betastar
			mu[i]=mustar
		} else {
			beta[i]=beta[i-1]
			mu[i]=mu[i-1]
		}

		#Multiple-try MH step on sigma2 using
		#normal random walk on U=log(sigma2)
		w1=rep(0,r)
		w2=rep(0,r)
		Usig=log(sigma2[i-1])
		Usigstar=rnorm(r,Usig,omega_sigma)
		sigma2star=exp(Usigstar)
		for (j in 1:r)	{
    			w1[j]=prod(1/sqrt(2*pi*sigma2star[j])*exp(-(theta[i,]-beta[i])^2/(2*sigma2star[j])))*dgamma(sigma2star[j],A,B)*exp(Usigstar[j])*dnorm(Usig,Usigstar[j],omega_sigma)*(dnorm(Usig,Usigstar[j],omega_sigma)*dnorm(Usigstar[j],Usig,omega_sigma))^(-1)
		}
		w1[is.na(w1)]=0
		if (sum(w1)!=0)	{
   			w=w1/sum(w1)
    			Usigstarj=sample(Usigstar,1,prob=w)
    			Usigdstar=rnorm((r-1),Usigstar,omega_sigma)
    			Usigdstar=c(Usigdstar,Usig)
    			sigma2dstar=exp(Usigdstar)
    			for (j in 1:r)	{
        			w2[j]=prod(1/sqrt(2*pi*sigma2dstar[j])*exp(-(theta[i,]-beta[i])^2/(2*sigma2dstar[j])))*dgamma(sigma2dstar[j],A,B)*exp(Usigdstar[j])*dnorm(Usigstarj,Usigdstar[j],omega_sigma)*(dnorm(Usigdstar[j],Usigstarj,omega_sigma)*dnorm(Usigstarj,Usigdstar[j],omega_sigma))^(-1)
        		}
    			w2[is.na(w2)]=0
    			if (sum(w2)!=0) {R=sum(w1)/sum(w2)} else {R=0}
		} else {R=0}

		sigma2[i]=ifelse(runif(1)<min(1,R),exp(Usigstarj),sigma2[i-1])

		#M-H step on U using uniform random walk
		Ustar=U[i-1]+sample(seq(-omega_U,omega_U,1),1)
		if (Ustar*k<T)	{R=0} else {
			fUstar<-dbinom(T,k*(Ustar),mu[i])*1/Ustar
			fU<-dbinom(T,k*(U[i-1]),mu[i])*1/U[i-1]
			R<-fUstar/fU
		}
		U[i]<-ifelse(runif(1) < min(1,R),Ustar,U[i-1])
		N[i]=U[i]+n[i]



	} #end chain loop

	out=list(N,theta,beta,sigma2,psi,mu,n,Q)
	return(out)
}