Skip to content

density function always evaluating to zero

4 messages · David Winsemius, napps22, Sarah Goslee

napps22
# 
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?

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.
David Winsemius
#
On Dec 3, 2011, at 5:42 PM, napps22 wrote:

            

      
      
      
    

        
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.

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      napps22
    

    

      
      #
    

  


  

      
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.

  
  


  


      

    

    

      
      

        


  

        

      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: