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.app...@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
______________________________________________
R-help@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.
