density function always evaluating to zero
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