归一化积分函数中的错误:“length = 21 in coercion to logical(1)”

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

提问人:david clarck 提问时间:9/16/2023 最后编辑:mkrieger1david clarck 更新时间:9/17/2023 访问量:33

问:

我在数值积分方面遇到了一个问题,无法对后验分布进行归一化。我的代码是:

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

出现的错误是

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

我该如何解决这个问题?

r 值积分 对数似然

评论


答:

1赞 jblood94 9/16/2023 #1

prior需要矢量化:

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

但请注意,由于在 [, ] 上是均匀的,并且积分的范围也是 [, ],只是 ,并且归一化常数是 。此外,由于积分的范围以均值为中心,因此积分最终是似然 () 中正态分布的均值:priora1 - sigma_a1a1 + sigma_a1a1 - sigma_a1a1 + sigma_a1posteriordnorm(a_true, a1, sigma_a1)/(2*sigma_a1)diff(pnorm(a1 + c(1, -1)*sigma_a1, a1, sigma_a1))/(2*sigma_a1)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

换言之,所有这些似乎都在计算均值和标准差的正态分布的均值,在 和 处截断。但是对称地截断尾巴不会改变平均值。a1sigma_a1a1 - sigma_a1a1 + sigma_a1

评论

0赞 david clarck 9/16/2023
非常感谢你的解释,太棒了。