density function always evaluating to zero

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

density function always evaluating to zero

Thread (4 messages)