# [R] compute posterior mean by numerical integration

Andreas Wittmann andreas_wittmann at gmx.de
Sat Sep 27 12:30:19 CEST 2008

```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

```

More information about the R-help mailing list