[R] density function always evaluating to zero
David Winsemius
dwinsemius at comcast.net
Sun Dec 4 05:50:15 CET 2011
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
More information about the R-help
mailing list