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.
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 at 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
Sorry that was my poor copying and pasting. Here's the correct R code. The
problem does seem to be with the function I define as f.
# Model selection example in a bayesian framework
# two competiting non-nested models
# M0: y_t = alpha * x1^2 + e_t
# M1: y_t = beta * x1^4 + e_t
# where e_t ~ iidN(0,1)
#generate data
x1 <- runif(100, min = -10, max = 10)
y <- 2 * x1^2 + rnorm(100)
n <- length(y)
k <- 1
v <- n - k
# want the posterior probabilities of the two models
# postprob_mj = p(y|model j true) * priorprob_mj / p(y)
# need to calculate the integral of p(y|alpha,h)p(alpha,h)
# and the integral of p(y|beta,h)p(beta,h)
# recall we chose p(alpha,h) = p(beta,h) = 1/h
# need to sample from the posterior to get an approximation of the integral
# # # # # # # # 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
# 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 m0
f <- function(alpha,h) {
e <- y - alpha * x1^2
const <- (2*pi)^(-n/2) / ( gamma(v/2) * (v*s2_m0/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-tp4155181p4156746.html
Sent from the R help mailing list archive at Nabble.com.