rnumerical-integrationlog-likelihood

Error in normalized integrate function: "length = 21 in coercion to logical(1)"


I have a problem with a numerical integration to normalize a posteriori distribution. My code is:

a1 <- 0.012
sigma_a1 <- 0.001


likelihood <- function(a_true) {
  (1 / (sqrt(2 * pi) * sigma_a1 )) * 
    exp(-((a1- a_true) ^ 2) / (2 * (sigma_a1 ^ 2)))
}

prior <- function(a_true) {
  if (a_true >= (a1 - sigma_a1 ) && a_true <= (a1+ sigma_a1 )) {
    return (1 / (2 * sigma_a1 ))
  } else {
    return (0)
  }
}

posterior <- function(a_true) {
  likelihood(a_true) * prior(a_true)
}

# step: Normalization of Posterior Distribution
integr <- integrate(posterior, lower = a1 - sigma_a1 , upper = a1 + sigma_a1)

a_true_estimated <- integrate(function(a_true) {
  a_true* posterior(a_true) / integral$value
}, lower = a1 - sigma_a1 , upper = a1 - sigma_a1)$value

The error presented was

Error in z_verdadeiro >= (z_photo - sigma_z_galphoto) && z_verdadeiro <= :
'length = 21' in coercion to 'logical(1)'

How can I fix this?


Solution

  • prior needs to be vectorized:

    a1 <- 0.012
    sigma_a1 <- 0.001
    
    
    likelihood <- function(a_true) {
      (1 / (sqrt(2 * pi) * sigma_a1 )) * 
        exp(-((a1- a_true) ^ 2) / (2 * (sigma_a1 ^ 2)))
    }
    
    prior <- function(a_true) {
      (a_true >= (a1 - sigma_a1) & a_true <= (a1 + sigma_a1))/(2*sigma_a1)
    }
    
    posterior <- function(a_true) {
      likelihood(a_true) * prior(a_true)
    }
    
    # step: Normalization of Posterior Distribution
    integral <- integrate(posterior, lower = a1 - sigma_a1 , upper = a1 + sigma_a1)
    
    a_true_estimated <- integrate(function(a_true) {
      a_true* posterior(a_true) / integral$value
    }, lower = a1 - sigma_a1 , upper = a1 + sigma_a1)$value
    
    a_true_estimated
    #> [1] 0.012
    

    But notice that since prior is uniform over [a1 - sigma_a1, a1 + sigma_a1], and the extents of the integration are also [a1 - sigma_a1, a1 + sigma_a1], posterior is simply dnorm(a_true, a1, sigma_a1)/(2*sigma_a1), and the normalization constant is diff(pnorm(a1 + c(1, -1)*sigma_a1, a1, sigma_a1))/(2*sigma_a1). Furthermore, since the extents of the integration are centered on the mean, the integral ends up being the mean of the normal distribution in the likelihood (a1):

    a1 <- 0.012
    sigma_a1 <- 0.001
    
    a_true_estimated <- integrate(
      function(a_true) a_true*dnorm(a_true, a1, sigma_a1),
      lower = a1 - sigma_a1 , upper = a1 + sigma_a1
    )$value/diff(pnorm(a1 + c(-1, 1)*sigma_a1, a1, sigma_a1))
    
    a1 - a_true_estimated
    #> [1] 1.734723e-18
    

    In other words, all this appears to be calculating the mean of a normal distribution with mean a1 and standard deviation sigma_a1, truncated at a1 - sigma_a1 and a1 + sigma_a1. But truncating the tails symmetrically does not change the mean.