#####################################################################################
# Code for: Jeff Gill and George Casella. "Nonparametric Priors For Ordinal         #
# Bayesian Social Science Models: Specification and Estimation. Journal of the      #
# American Statistical Association, Forthcoming JASA, 2009.                         #
#                     		       Copyright 2008, Jeff Gill and George Casella #
#####################################################################################

# PRE-LOOP EFFICIENCIES
tau      <- rgamma(1,shape=a,scale=b)
mu       <- rnorm(1,0,d*tau^2)
sig.tau  <- sqrt(sigma^2 + tau^2)
A.inv    <- solve(A)
sigma.sq <- sigma^2
tau.sq   <- tau^2
beta     <- rbind( beta,matrix(rep(NA,num.reps*ncol(beta)),ncol=ncol(beta)) )
theta    <- rbind( theta,matrix(rep(theta,num.reps),byrow=TRUE,ncol=ncol(theta)) )
U        <- matrix(NA,nrow=num.reps,ncol=n)

# MCMC LOOPING
for (M in 1:num.reps)  {
    if (M %% 10 == 0) print(M)

    #print("DRAW u, RTNORM")
    for (i in 1:n)   	
	u[i] <- rtrun( mu=X[i,]%*%beta[M,]+psi[i], sigma=sigma.sq, a=theta[M,which.max(Y[i,])], b=theta[M,(which.max(Y[i,])+1)] )

    #print("GET C.vec AND UPDATE C.mat WITH MORE ASSOCIATIONS")
    for (i in 1:n)  {
	nil[i]     <- length(C.vec[C.vec==C.vec[i]])-1
	f.y.psi[i] <- prod( (pnorm((theta[M,-1]     - X[i,]%*%beta[M,] - psi[i])/sigma) - 
                             pnorm((theta[M,-C.num] - X[i,]%*%beta[M,] - psi[i])/sigma))^Y[i,] )
	H[i]       <- pnorm( (theta[M,(which.max(Y[i,])+1)] - mu + X[i,]%*%beta[M,])/sig.tau ) -
	              pnorm( (theta[M,(which.max(Y[i,])+0)] - mu + X[i,]%*%beta[M,])/sig.tau ) 
	b          <- sum(nil[i]*f.y.psi[i])/(n-1+m) + length(nil[nil==0])*m/(n-1+m)
	C.prob     <- nil*f.y.psi/(n-1+m);  C.prob[C.prob==0] <- m*H[i]/(n-1+m); C.prob <- C.prob/sum(C.prob)
	C.vec[i]   <- sample(x=1:n,size=1,replace=FALSE,prob=C.prob)
	if (duplicated(c(C.vec[-i],C.vec[i]))[n] == FALSE)  {
		T      <- u[i]-X[i,]%*%beta[M,]
		psi[i] <- rnorm(1, (T*tau.sq + mu*sigma.sq)/sig.tau, ((sigma.sq)*(tau.sq))/sig.tau )
	}
    }
    D.mat <- matrix(0,nrow=n,ncol=n)
    D.mat[cbind(C.vec,1:n)] <- 1
    D.mat <- t(D.mat) %*% D.mat
    D.mat[row(D.mat) <= col(D.mat)] <- 0
    C.mat <- C.mat + D.mat
    C.summary[M,1] <- length(unique(C.vec)); C.summary[M,2] <- mean(table(C.vec))      

    #print("UPDATE psi VALUES, USES test.vec AS A SCREEN")
    done.vec <- rep(0,n)
    for (i in 1:n)  {
	if (done.vec[i] == 0)  {
	    test.vec <- C.vec; test.vec[test.vec != C.vec[i]] <- 0; test.vec[test.vec == C.vec[i]] <- 1
	    done.vec <- done.vec + test.vec
	    T        <- sum(test.vec*(u-X%*%beta[M,]))/sum(test.vec)
	    psi      <- psi*(1-test.vec) +  test.vec*rnorm(1, (sum(test.vec)*T*tau.sq + mu*sigma.sq)/(sigma.sq + sum(test.vec)*tau.sq), 
							      ((sigma.sq)*(tau.sq))/(sigma.sq + sum(test.vec)*tau.sq))
	}
    }

    #print("SAMPLE betaS AND thetaS")
    B            <- ( sigma.beta.sq*t(X)%*%(u-psi) + mu.beta*sigma.sq )/( sigma.sq * sigma.beta.sq )
    beta[(M+1),] <- rmultinorm(1,A.inv%*%B,A.inv) 
    bounds.vec   <- Y*u;  bounds.vec[bounds.vec==0] <- NA;  bounds.vec <- apply(bounds.vec,2,range,na.rm=TRUE)
    for (i in 1:(ncol(bounds.vec)-1))  theta[(M+1),(i+1)] <- runif(1,min=bounds.vec[2,i],max=bounds.vec[1,(i+1)])

    #print("UPDATE mu AND tau, c=0 FOR NOW")
    r       <- length(table(psi))
    tau     <- 1/rgamma(1, shape=(r+1)/2, scale=sum((table(psi)-mu)^2)/2 - ((mu-0)^2)/(2*d) + 1/b)
    tau.sq  <- tau^2
    sig.tau <- sqrt(sigma.sq + tau.sq)
    mu      <- rnorm(1, ((d*tau.sq)/(r*d*tau.sq+1))*sum(table(psi))+e/d, (d*tau.sq)/(r*d*tau.sq+1) )

    #print("UPDATE m PARAMETER, ASSUMED INTEGER WITH LOG-UNIFORM PRIOR, APPROXIMATION")
    c.star <- length(unique(C.vec))
    for (i in 1:length(m.val.vec))  {
	log.sum <- sum(log(n:i))
	m.prob.vec[i] <- c.star * log(m.val.vec[i]) - log.sum
    }
    m.prob.vec <- m.prob.vec - sum(m.prob.vec)
    m.prob.vec <- m.prob.vec/sum(m.prob.vec)
    m <- sample(m.val.vec,size=1,prob=m.prob.vec)
}