################################################################################################## # IMPORTANCE SAMPLING ROUTINES FOR S/R # This is a _very_ simple importance sampling routine mainly for use by my students, # but others are welcome to use/modify/etc. provided that it is understood that there # exists no warranty or guarantee of any sort. # # Draw a random sample from the multivariate normal or t pdf with mean equal # to the mle estimates from the input log-likelihood function, and the standard error # equal to the sqrt of the diagonal of the generalized chol/inv var/covar matrix. # The function "sir" performs importance sampling with a normal approximation # as a default, but one can easily change the specification by stating a df. # Actually its important to note that this really implements Rubin's SIR since # that way I have an empirical maximum to divide by. Please email if you find # heinous errors or the like. -Jeff Gill. ################################################################################################## # Call: # sir(data.mat, the original data matrix organized in columns # theta.vector, the initial coefficient estimates # theta.matrix, the initial vc matrix # M, the number of draws # m, the desired number of accepted values # tol, the rounding/truncing tolerance # ll.func, the likelihood function for the empirical posterior. # df) the df for using the t distribution as the approximation distribution, # default is 0 for gaussian. ################################################################################################## sir <- function(data.mat,theta.vector,theta.mat,M,m,tol=1e-06,ll.func,df=0) { importance.ratio <- rep(NA,M) rand.draw <- rmultnorm(M,theta.vector,theta.mat,tol = 1e-04) if (df > 0) rand.draw <- rand.draw/(sqrt(rchisq(M,df)/df)) empirical.draw.vector <- apply(rand.draw,1,ll.func,data.mat) if (sum(is.na(empirical.draw.vector)) == 0) { print("SIR: finished generating from posterior density function") print(summary(empirical.draw.vector)) } else { print(paste("SIR: found",sum(is.na(empirical.draw.vector)), "NA(s) in generating from posterior density function, quiting")) return() } if (df == 0) { normal.draw.vector <- apply(rand.draw,1,normal.posterior.ll,data.mat) } else { theta.mat <- ((df-2)/(df))*theta.mat normal.draw.vector <- apply(rand.draw,1,t.posterior.ll,data.mat,df) } if (sum(is.na(normal.draw.vector)) == 0) { print("SIR: finished generating from approximation distribution") print(summary(normal.draw.vector)) } else { print(paste("SIR: found",sum(is.na(normal.draw.vector)), "NA(s) in generating from approximation distribution, quiting")) return() } importance.ratio <- exp(empirical.draw.vector - normal.draw.vector) importance.ratio[is.finite=F] <- 0 importance.ratio <- importance.ratio/max(importance.ratio) if (sum(is.na(importance.ratio)) == 0) { print("SIR: finished calculating importance weights") print(summary(importance.ratio)) } else { print(paste("SIR: found",sum(is.na(importance.ratio)), "NA(s) in calculating importance weights, quiting")) return() } accepted.mat <- rand.draw[1:2,] while(nrow(accepted.mat) < m+2) { rand.unif <- runif(length(importance.ratio)) accepted.loc <- seq(along=importance.ratio)[(rand.unif-tol) <= importance.ratio] rejected.loc <- seq(along=importance.ratio)[(rand.unif-tol) > importance.ratio] accepted.mat <- rbind(accepted.mat,rand.draw[accepted.loc,]) rand.draw <- rand.draw[rejected.loc,] importance.ratio <- importance.ratio[rejected.loc] print(paste("SIR: cycle complete,",(nrow(accepted.mat)-2),"now accepted")) } accepted.mat[3:nrow(accepted.mat),] } ##################################################################################################