[R] density function always evaluating to zero

[David Winsemius]

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