[R] compute posterior mean by numerical integration

[Andreas Wittmann]

Dear R useRs,

i try to compute the posterior mean for the parameters omega and beta 
for the following
posterior density. I have simulated data where i know that the true 
values of omega=12
and beta=0.01. With the function postMeanOmega and postMeanBeta i wanted 
to compute
the mean values of omega and beta by numerical integration, but instead 
of omega=12
and beta=0.01 i get omega=11.49574 and beta=1.330105. I don't know what 
is wrong
in my implementation, i guess the computation by numerical integration 
strongly depends
on the integration limits low, upw, lob and upb, but what are good 
choices for these to get
reasonable results?


#######################################################

## simulated data with omega=12, beta=0.01
data <- c(8, 1, 6, 14, 1, 0, 16, 37, 15, 17, 6, 57)
 
## log likelihood function
"loglike" <- function(t, omega, beta)
{
  n <- length(t)-1
  res <- n * log(omega) + n * log(beta) - beta * 
sum(cumsum(t[-length(t)])) -
         omega * (1-exp(-beta * cumsum(t)[n]))
         
  return(res)
}
 
## log prior density
"prior" <- function(omega, beta, o1, o2, b1, b2)
{
  if(o1 && o2 && b1 && b2)
    res <- dgamma(omega, o1, o2, log=T) + dgamma(beta, b1, b2, log=T)
  else
    res <- 0 ## noninformative prior
    
  return(res)
}
 
## log posterior density
"logPost" <- function(t, omega, beta, o1, o2, b1, b2)
{
  res <- loglike(t, omega, beta) + prior(omega, beta, o1, o2, b1, b2)
   
  return(res)
}
 
## posterior normalizing constant
"PostNormConst" <- function(t, low, upw, lob, upb, o1, o2, b1, b2)
{
  "g" <- function(beta, ...)
  {
    integrate(function(omega) {logPost(t=t, omega, beta, o1=o1, o2=o2, 
                                       b1=b1, b2=b2)}, low, upw)$value
  }
 
  res <- integrate(Vectorize(function(beta) {g(beta, low=low, upw=upw, t=t,
                                               o1=o1, o2=o2, b1=b1, 
b2=b2)}), 
                                               lob, upb)$value
 
  return(res)
}
 
## posterior mean omega
"postMeanOmega" <- function(t, norm, low, upw, lob, upb, o1, o2, b1, b2)
{
  "g" <- function(beta, ...)
  {
    integrate(function(omega) {logPost(t=t, omega, beta, o1=o1, o2=o2,
                                       b1=b1, b2=b2) * omega}, low, 
upw)$value
  }
 
  tmp <- integrate(Vectorize(function(beta) {g(beta, low=low, upw=upw, t=t,
                                               o1=o1, o2=o2, b1=b1, 
b2=b2)}), 
                                               lob, upb)$value
 
  res <- tmp / norm
 
  return(res)
}
 
## posterior mean beta
"postMeanBeta" <- function(t, norm, low, upw, lob, upb, o1, o2, b1, b2)
{
  "g" <- function(beta, ...)
  {
    integrate(function(omega) {logPost(t=t, omega, beta, o1=o1, o2=o2,
                                       b1=b1, b2=b2) * beta}, low, 
upw)$value
  }
 
  tmp <- integrate(Vectorize(function(beta) {g(beta, low=low, upw=upw, t=t,
                                               o1=o1, o2=o2, b1=b1, 
b2=b2)}), 
                                               lob, upb)$value
 
  res <- tmp / norm
 
  return(res)
}
 
low <- 3; upw <- 20;
lob <- 0; upb <- 2;
 
norm <- PostNormConst(t=data, low=low, upw=upw, lob=lob, upb=upb, o1=0, 
o2=0, b1=0, b2=0)
postMeanOmega(t=data, norm=norm, low=low, upw=upw, lob=lob, upb=upb, 
o1=0, o2=0, b1=0, b2=0)
postMeanBeta(t=data, norm=norm, low=low, upw=upw, lob=lob, upb=upb, 
o1=0, o2=0, b1=0, b2=0)

#######################################################


If you have any advice, i would be very thankful,

best regards

Andreas