multivspost2blogandtausqtot: Function to evaluate marginalized and reparameterized...

In CARrampsOcl: Reparameterized and marginalized posterior sampling for conditional autoregressive models, OpenCL implementation

Description

Function to evaluate marginalized and reparameterized posterior.

Usage

multivspost2blogandtausqtot(smat, alpha, beta, D, y, By, k)

Arguments

smat
alpha
beta
D
y
By
k

Value

Matrix in which the first column is the log of the marginalized posterior density of each set of candidate values and the second column is a draw from the posterior conditional density of tausqtot given the corresponding set of candidate values.

Author(s)

Kate Cowles

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function(smat, alpha, beta, D, y, By, k )

{
# function to evaluate log marginal posterior of s and generate tausqtot

smat <- as.matrix(smat)

n <- length(y)

F <- ncol(D) + 1

logpostvect <- numeric()
tausqtot <- numeric()

sumalpha <- sum(alpha)


# change to use apply 09/29/09 MKC

genst <- function( s )
{
        s0 <- 1-sum(s)

        neweigennumer <- s[1] * D[,1]
        if( F > 2)
        for(j in 2:(F-1))
                neweigennumer <- neweigennumer + s[j] * D[,j]

        neweigendenom <-  neweigennumer + s0

        neweigen <- s0 * neweigennumer / neweigendenom
# corrected to (alpha - 1)  09/18/09
        #logpostdensnumer <- sum( log(c(s0,s)) * alpha  ) + 
        logpostdensnumer <- sum( log(c(s0,s)) * (alpha-1)  ) + 
                        sum( log(neweigen[ neweigen > 0 ]) ) / 2
        whole <- sum( neweigen * By^2)

        newbeta <- whole/2 +  sum( c(s0,s) * beta )
        newalpha <- (sumalpha + (n-k)/2 )
        logpostdensdenom <-  log( newbeta)* newalpha
        #logpostvect[i] <- logpostdensnumer - logpostdensdenom
        #tausqtot[i] <- rgamma(1, newalpha, newbeta)
        c( logpostdensnumer - logpostdensdenom, 
        rgamma(1, newalpha, newbeta))
  }


if( multicoreflag )
    output <- t( parApply(cl, smat, 1, genst ) )
else
    output <- t( apply( smat, 1, genst ) )
output
 
  }

CARrampsOcl documentation built on May 2, 2019, 3:27 a.m.