[R] density function always evaluating to zero

[Sarah Goslee]

I didn't try out your extensive code, but here's one potentially serious
problem:

You only pass two arguments to f(), alpha and h, but within f you nonetheless
use x1 and y and several other things. This is bad practice, and dangerous:
you should pass all the necessary arguments to f(), not rely on the environment
to contain them.

After you fix that, working through your function step by step at the command
line and looking at the results of each line of code can be very educational.

Sarah

On Sun, Dec 4, 2011 at 7:47 AM, napps22 <n.j.apps22 at gmail.com> wrote:
> 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
>
>
>

-- 
Sarah Goslee
http://www.functionaldiversity.org