Dear R users,
I'm trying to carry out monte carlo integration of a posterior density
function which is the product of a normal and a gamma distribution. The
problem I have is that the density function always returns 0. How can I
solve this problem?
Here is my code
#generate data
x1 <- runif(100, min = -10, max = 10)
y <- 2 * x1^2 + rnorm(100)
# # # # # # # # Model 0 # # # # # # #
z <- x1^2
M <- sum(z^2)
MI <- 1/M
zy <- crossprod(z,y)
alpha.ols <- MI * zy
resid_m0 <- y - z * alpha.ols
s2_m0 <- sum(resid_m0^2)/v
# use gibbs sampler to sample from posterior densities
n <- length(y)
k <- 1
v <- n - k
# set up gibbs sampler
nrDraws <- 10000
h_sample_m0 <- rgamma(nrDraws, v/2, v*s2_m0/2)
alpha_sample <- matrix(0, nrow = nrDraws, ncol = 1)
for(i in 1:nrDraws) {
alpha_sample[i] <- rnorm(1,alpha.ols,(1/h_sample_m0[i]) * MI)
}
# define posterior density for model 0
f <- function(alpha,h) {
e <- y - alpha * x1^2
const <- (2*pi)^(-n/2) / ( gamma(v/2) * (v*s2_m1/2)^(-v/2) )
kernel <- h^((v-1)/2) * exp((-(h/2) * sum(e^2)) - (v*s2_m0)*h)
x <- const * kernel
return(x)
}
# calculate approximation of p(y|m_0)
m0 <-f(alpha_sample,h_sample_m0)
post_m0 <- sum(m0) / nrDraws
--
View this message in context:
http://r.789695.n4.nabble.com/density-function-always-evaluating-to-zero-tp4155181p4155181.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
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.
