R/probitFit.r

In dhelkey/dgrank: Construct Draper-Gittoes Rankings of Institutions

probitFit = function(y, X,  prior_var_vec, iters = 100){
	#' MCMC Samples for Probit Regression
	#'
	#' \code{probitFit} generates MCMCM samples of Probit regression, following
	#' the data-augmentation Gibbs sampler of Albert and Chib (1993).
	#' Rquires truncnorm package to sample from truncated normal.
	#' Assumes prior centered at 0, and diagonal prior covariance structure.
	#' Assumes response y is coded as 0-1.
	#'
	#' @param y <nX1> observed data vector
	#' @param X <nXp> design matrix (takes objects of class matrix or sparce Matrix)
	#' @param prior_var_vec <px1> prior variance vector
	#' @param iters Integer; number of desired output iterations

	#Number of coefficients
	n = dim(X)[1]
	p = dim(X)[2]
	B_inv = diag(1/prior_var_vec, ncol = p) #inverse prior variance mat
	
	#Truncation limits, 1 row per observation, colums are lower and upper limits.
	trunc_limits = matrix(0, nrow = n, ncol = 2) 
	trunc_limits[y==0 ,1] = -Inf
	trunc_limits[y==1, 2] = Inf
		

	
	#Precompute Cholesky decomposition (to efficiently sample from multivariate normal)
	XtX = Matrix::crossprod(X)
	Xty = Matrix::crossprod(X,y)
	A = XtX + B_inv
	R = chol(A) #Cholesky decomposition: A = t(R) %*% R

	#Start MCMC algorithm at Least-Squares estimate (guarenteed to exist w/ nonzero prior variance)
	beta_ls = linCholSolver(R, Xty)

	#Generate the normal rancom variates outside of loop
	rnorm_mat = matrix( rnorm(p * iters, mean = 0, sd = 1), 
				ncol = p)
				
	#Helper function to sample Z|... from truncated normal
	sampleZ = function(beta_vec){
		#Samples z given beta from truncated normal distribtion
		truncnorm::rtruncnorm(n, trunc_limits[ ,1], trunc_limits[ ,2], 
							mean =  as.numeric(X %*% beta_vec))							
	}
	
	#MCMC storage, each row is an iteration
	beta_mat = matrix(0, nrow = iters, ncol = p )
	beta_mat[1, ] = as.numeric(beta_ls)
	z_latent = sampleZ(beta_mat[1, ])
	for (i in 2:iters){
		#Gibbs sample beta given z
		Xtz = Matrix::crossprod(X, z_latent)
		beta_z = linCholSolver(R,  Xtz)
		beta_mat[i, ] =  as.numeric(beta_z +
			backsolve(R, rnorm_mat[i, ]))

		#Gibbs sample z
		z_latent = sampleZ(beta_mat[i, ])
	}
	return(beta_mat)
}