# R/getLogMargLikEstimate.R In glmBfp: Bayesian Fractional Polynomials for GLMs

#### Documented in getLogMargLikEstimate

```#####################################################################################
## Author: Daniel Sabanés Bové [daniel *.* sabanesbove *a*t* ifspm *.* uzh *.* ch]
## Project: BFPs for GLMs
##
## Time-stamp: <[getLogMargLikEstimate.R] by DSB Mon 05/07/2010 13:36 (CEST)>
##
## Description:
## Compute the Chib-Jeliazkov log marginal likelihood estimate from the MCMC output
##
## History:
## 18/02/2010   separate code file;
##              correct bug for computing Omega(0): head(index, -0L) always gives
##              a length zero vector instead of the whole "index"!!!
#####################################################################################

##' Compute the Chib-Jeliazkov log marginal likelihood estimate from the MCMC output
##'
##' @param numeratorTerms terms for the sum in the numerator of the posterior density ordinate
##' estimate at the high density point
##' @param denominatorTerms terms for the sum in the denominator of the posterior density ordinate
##' estimate at the high density point
##' @param highDensityPointLogUnPosterior log unnormalized posterior (i.e. likelihood times prior)
##' of the high density point
##' @param endLag up to which lag should the covariance incorporate the uncertainty? (default: 40)
##' @return list with resulting \dQuote{estimate} and \dQuote{standardError}.
##'
##' @keywords internal
##'
getLogMargLikEstimate <- function(numeratorTerms,
denominatorTerms,
highDensityPointLogUnPosterior,
endLag=40L)
{
## note that here it is assumed that there are equally many numerator and
## denominator numbers, therefore
stopifnot(identical(length(numeratorTerms),
length(denominatorTerms)),
is.double(numeratorTerms),
is.double(denominatorTerms),
is.double(highDensityPointLogUnPosterior))

## first combine each numerator term with corresponding denominator term
h <- rbind(numeratorTerms,
denominatorTerms)

## the estimate for that is
hHat <- rowMeans(h)

## so the estimate for the log posterior density ordinate is:
logPosteriorDensityEstimate <- log(hHat[1]) - log(hHat[2])

## then compute the standard error for this estimate:
hDiffHat <- h - hHat

B <- ncol(hDiffHat)
index <- seq_len(B)

## function to compute the autocovariance estimate for a given lag
getOmega <- function(lag)
{
## get index for head part of time series
if(identical(lag, 0L))
index
else

## and for tail
tailIndex <-
if(identical(lag, 0L))
index
else
tail(index, - lag)

## now build up result
ret <- matrix(data=0,
nrow=2L,
ncol=2L)

for(l in seq_len(B - lag))
{
ret <- ret + tcrossprod(hDiffHat[, tailIndex[l]],
}

return(ret)
}

covMat <- getOmega(0L) / B

for(s in seq_len(endLag))
{
thisOmega <- getOmega(s)
covMat <- covMat + (1 - s / (endLag + 1)) * (thisOmega + t(thisOmega)) / B
}

## now covMat is ready to get delta ruled:
aHat <- c(1 / hHat[1], - 1 / hHat[2])
se <- sqrt(t(aHat) %*% covMat %*% aHat)

## return a list with the log marg lik estimate and the standard error
return(list(estimate=highDensityPointLogUnPosterior - logPosteriorDensityEstimate,
standardError=se))
}
```

## Try the glmBfp package in your browser

Any scripts or data that you put into this service are public.

glmBfp documentation built on Sept. 7, 2018, 3:01 a.m.