该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
(见这里),但它涉及赖特的广义合流超几何函数,它还没有在 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
问题:
- 你能做得比这更好吗?
- 是否有一些关于数值积分例程的东西,这些事情的专家可以对这里发生的事情有所了解?我有一个偷偷摸摸的怀疑,
t
余弦会迅速波动,这会导致问题......? - 是否存在更适合此类问题的现有 R 例程/程序包,有人可以指点我到一个合适的位置(在山上)开始攀登吗?
评论:
t
是的,用作函数参数是不好的做法。- 我计算了
shape > 1
使用 Maple 发布的结果的确切答案,以及被brute-force-integrate-by-the-definition-with-R
踢的 Maple 的屁股。也就是说,我可以在几分之一秒内得到相同的答案(达到数值精度),而价格甚至更小。
编辑:
我打算写下我正在寻找的确切积分,但似乎这个特定站点不支持 MathJAX,所以我将提供链接。我正在寻找数值评估Weibull 分布的特征函数以获得合理的输入(无论这意味着什么)。该值是一个复数,但我们可以将其拆分为实部和虚部,这就是我所说的及以上。t
Rp
Ip
最后一条评论:维基百科列出了 Weibull cf 的公式(无限级数),并且该公式与我在上面引用的论文中证明的公式相匹配,但是,该级数仅被证明适用于shape > 1
. 该案0 < shape < 1
仍然是一个悬而未决的问题;有关详细信息,请参阅论文。