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
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
