# [R] density function always evaluating to zero

napps22 n.j.apps22 at gmail.com
Sun Dec 4 13:47:45 CET 2011

```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

--
View this message in context: http://r.789695.n4.nabble.com/density-function-always-evaluating-to-zero-tp4155181p4156746.html
Sent from the R help mailing list archive at Nabble.com.

```