On Dec 3, 2011, at 5:42 PM, napps22 wrote:
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?
Your code throws errors due to several missing objects due in part to
a missing "v", but moving that higher in the code does not cure
everything. A missing 's2_m1' is also listed as a source of error.
You might want to try in a clean session and redo your debugging so
that you know what variables are in your workspace and their values.
The usual debugging method is to sprinkle print() statements in your
code to monitor where things might not be proceeding as planned.
--
David.
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.
David Winsemius, MD
West Hartford, CT
______________________________________________
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.
