在 R 中发布 abour 集成（）函数

Question

我有一个关于 R 中的积分（）函数的问题。当我使用积分函数求后验均值时，结果非常小（小于绝对误差）。如果我将上限从 Inf 更改为 0.2 或更小的数字，则结果数字更有意义。为什么较大的上限积分功能无法正常工作？

非常感谢！

以下是我的代码：

dat = read.table("c://AAPL_logret.txt")$V1
mu = 0.005
theta = 200

posterior = function (sigma2, dat)
{
  numerator = function (sigma2, dat)
  {
    out = dexp(sigma2,theta)
    for (i in 1:length(dat))
    {
     out = out * dnorm(dat[i],mu, sqrt(sigma2))
    }
    return(out)   

  }

  denominator = integrate(numerator,lower = 0, upper = Inf, dat=dat)$value

  return (numerator(sigma2,dat)/denominator)
}

curve(posterior(x,dat), from =0, to = 0.006,xlab=expression(sigma^2), 
      ylab = "Density", cex.axis = 1.3, cex.lab=1.3, col=4,lwd=1.5,lty=4,n=20000)

posterior_times_sigma2 = function(sigma2,dat)
{
  return(sigma2 * posterior(sigma2,dat))
}

posteriormean = integrate(posterior_times_sigma2, lower =0, upper= Inf, dat=dat)$value

Answer 1

您还没有发布“AAPL_logret.txt”文件的内容，所以我生成了可能看起来像这样的虚假数据，例如：

dat <- rlnorm(200, -4, 1)

有了这些数据，您的 curve() 图显示后验分布的大部分概率质量包含在 0 到 0.006 之间：

 1-integrate(posterior, lower=0, upper=0.006, dat=dat)$value

大约为 4e-11，因此后验在 0.006 之后与零在数值上无法区分。集成（）的 R 帮助页面解释了这里出了什么问题：

“像所有数值积分例程一样，这些函数评估 在有限的点集上起作用。如果函数是 在几乎所有的 范围，结果和误差估计可能是 严重错误。”

我相信这里没有通用的解决方案：在任何给定的 您需要探索不同的选择。在你的情况下， 一种可能性是限制积分间隔：

integrate(posterior_times_sigma2, lower=0, upper=Inf, dat=dat)$value
integrate(posterior_times_sigma2, lower=0, upper=0.2, dat=dat)$value
integrate(posterior_times_sigma2, lower=0, upper=0.1, dat=dat)$value
integrate(posterior_times_sigma2, lower=0, upper=0.006, dat=dat)$value

不同的集成方法也可能有所帮助，例如 pracma 包提供的方法：

library(pracma)
integral(posterior_times_sigma2, xmin=0, xmax=Inf, dat=dat)

这是整个脚本：

##dat = read.table("c://AAPL_logret.txt")$V1

set.seed(1233)
dat <- rlnorm(100, -3, 0.1)

mu= 0.005
theta = 200

posterior = function (sigma2, dat){
numerator = function (sigma2, dat)  {
    out = dexp(sigma2, theta)
    for (i in 1:length(dat))    {
    out = out * dnorm(dat[i], mu, sqrt(sigma2))
    }
    return(out)   
}

denominator = integrate(numerator,lower = 0, upper = Inf, dat=dat)$value
return (numerator(sigma2,dat)/denominator)
}

curve(posterior(x,dat), from=0, to = 0.006,
  xlab=expression(sigma^2), ylab = "Density", cex.axis = 1.3,
  cex.lab=1.3, col=4, lwd=1.5, lty=4, n=20000)

integrate(posterior, lower=0, upper=0.006, dat=dat)$value

posterior_times_sigma2 = function(sigma2,dat){
return(sigma2 * posterior(sigma2,dat))
}

integrate(posterior_times_sigma2, lower=0, upper=Inf, dat=dat)$value
integrate(posterior_times_sigma2, lower=0, upper=0.2, dat=dat)$value
integrate(posterior_times_sigma2, lower=0, upper=0.1, dat=dat)$value
integrate(posterior_times_sigma2, lower=0, upper=0.006, dat=dat)$value



library(pracma)
integral(posterior_times_sigma2, xmin=0, xmax=Inf, dat=dat)
integral(posterior_times_sigma2, xmin=0, xmax=Inf, dat=dat, method="Kronrod")