数值整合一个棘手的功能

Question

该prob软件包以数字方式评估基本R分布的特征函数.几乎所有的发行版都有现有的公式.但是,在少数情况下,没有封闭形式的解决方案.例证:威布尔分布(但见下文).

对于Weibull特征函数,我基本上计算两个积分并将它们放在一起:

fr <- function(x) cos(t * x) * dweibull(x, shape, scale)
fi <- function(x) sin(t * x) * dweibull(x, shape, scale)
Rp <- integrate(fr, lower = 0, upper = Inf)$value
Ip <- integrate(fi, lower = 0, upper = Inf)$value
Rp + (0+1i) * Ip

是的,这很笨拙,但效果出奇的好!...... 啊哈,大部分的时间.用户最近报告了以下情况:

cfweibull(56, shape = 0.5, scale = 1)

Error in integrate(fr, lower = 0, upper = Inf) : 
  the integral is probably divergent

现在,我们知道积分不是发散的,所以它必然是一个数值问题.随着一些摆弄,我可以得到以下工作:

fr <- function(x) cos(56 * x) * dweibull(x, 0.5, 1)

integrate(fr, lower = 0.00001, upper = Inf, subdivisions=1e7)$value
[1] 0.08024055

没关系,但这不太对,加上它需要一些不能很好地扩展的摆弄.我一直在调查这个以获得更好的解决方案.我发现了最近发表的关于特征函数的"封闭形式" scale > 1(见这里),但它涉及Wright的广义汇合超几何函数,这在R(尚未)中没有实现.我查看档案中的integrate备选方案,那里有很多东西,看起来并不是很好.

作为搜索的一部分,我想通过反正切将积分区域转换为有限区间,瞧!看看这个:

cfweibull3 <- function (t, shape, scale = 1){
  if (shape <= 0 || scale <= 0) 
    stop("shape and scale must be positive")
  fr <- function(x) cos(t * tan(x)) * dweibull(tan(x), shape, scale)/(cos(x))^2
  fi <- function(x) sin(t * tan(x)) * dweibull(tan(x), shape, scale)/(cos(x))^2
  Rp <- integrate(fr, lower = 0, upper = pi/2, stop.on.error = FALSE)$value
  Ip <- integrate(fi, lower = 0, upper = pi/2, stop.on.error = FALSE)$value
  Rp + (0+1i) * Ip
}

> cfweibull3(56, shape=0.5, scale = 1)
[1] 0.08297194+0.07528834i

问题:

1. 你能做得比这更好吗？
2. 是否存在关于数字集成例程的事情,那些专注于此类事情的人可以对这里发生的事情有所了解？我有一种潜在的怀疑,即大量t的余弦会迅速波动而导致问题......？
3. 是否存在更适合此类问题的现有R例程/包,并且有人可以指向一个位置良好的位置(在山上)以开始攀爬吗？

评论:

1. 是的,使用t函数参数是不好的做法.
2. 我计算了shape > 1使用Maple发布的结果和brute-force-integrate-by-the-definition-with-R踢出的Maple屁股的确切答案.也就是说,我在一小部分时间内得到了相同的答案(达到数值精度),甚至是价格的一小部分.

编辑:

我打算记下我正在寻找的确切积分,但似乎这个特定网站不支持MathJAX所以我会给出链接.我正在寻找数值评估Weibull分布的特征函数以获得合理的输入(无论这意味着什么).这个值是一个复数,但是我们可以将它分成它的实部和虚部,这就是我所调用的和以上.tRpIp

最后一条评论:维基百科有一个公布的(无限系列)公式为Weibull cf,该公式与我在上面引用的论文中证明的公式匹配,但是,该系列只被证明适用shape > 1.案件0 < shape < 1仍然是一个悬而未决的问题; 请参阅论文了解详情.

Answer 1

您可能有兴趣看一下本文,它讨论了高振荡积分的不同积分方法 - 这就是您基本上想要计算的内容:http://citeseerx.ist.psu.edu/viewdoc/summary？ doi = 10.1. 1.8.6944

另外,另一个可能的建议是,你可能想要指定一个更小的,而不是无限制,因为如果你指定你想要的精度,那么根据weibull的cdf,你可以很容易地估计你可以有多少尾巴截短.如果你有一个固定的限制,那么你可以准确地(或几乎)指定细分的数量(例如,为了每个周期有几个(4-8)点).