sir {BaM}

sir

Description

Implementation of Rubin's SIR, see pages 338-341 (2nd Edition)

Usage

sir(data.mat,theta.vector,theta.mat,M,m,tol=1e-06,ll.func,df=0)

Arguments

data.mat

A matrix with two columns of normally distributed data

theta.vector

The initial coefficient estimates

theta.mat

The initial vc matrix

M

The number of draws

m

The desired number of accepted values

tol

The rounding/truncing tolerance

ll.func

loglike function for empirical posterior

df

The df for using the t distribution as the approx distribution

Author(s)

Jeff Gill

Examples

## Not run:  
sir <- function(data.mat,theta.vector,theta.mat,M,m,tol=1e-06,ll.func,df=0) {
 importance.ratio <- rep(NA,M)
 rand.draw <- rmultinorm(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),]
}
# The following are log likelihood functions that can be plugged into the sir function above.

logit.posterior.ll <- function(theta.vector,X) {
  Y <- X[,1]
  X[,1] <- rep(1,nrow(X))
  sum( -log(1+exp(-X
                  -log(1+exp(X)))))
}

normal.posterior.ll <- function(coef.vector,X) {
  dimnames(coef.vector) <- NULL
  Y <- X[,1]
  X[,1] <- rep(1,nrow(X))
  e <- Y - X
  sigma <- var(e)
  return(-nrow(X)*(1/2)*log(2*pi)
         -nrow(X)*(1/2)*log(sigma)
         -(1/(2*sigma))*(t(Y-X)*(Y-X)))
}

t.posterior.ll <- function(coef.vector,X,df) {
  Y <- X[,1]
  X[,1] <- rep(1,nrow(X))
  e <- Y - X
  sigma <- var(e)*(df-2)/(df)
  d <- length(coef.vector)
 return(log(gamma((df+d)/2)) - log(gamma(df/2))
       - (d/2)*log(df)
       -(d/2)*log(pi) - 0.5*(log(sigma))
       -((df+d)/2*sigma)*log(1+(1/df)*
                               (t(Y-X*(Y-X)))))
}

probit.posterior.ll <- function (theta.vector,X,tol = 1e-05) {
  Y <- X[,1]
  X[,1] <- rep(1,nrow(X))
  Xb <- X
  h <- pnorm(Xb)
  h[h<tol] <- tol
  g <- 1-pnorm(Xb)
  g[g<tol] <- tol
  sum( log(h)*Y + log(g)*(1-Y) )
}

## End(Not run)