13

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

问题:

  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仍然是一个悬而未决的问题;有关详细信息,请参阅论文。

4

4 回答 4

6

您可能有兴趣查看这篇论文,该论文讨论了高度振荡积分的不同积分方法——这就是您本质上要计算的内容: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1。 1.8.6944

另外,另一个可能的建议是,您可能希望指定一个较小的限制,而不是无限限制,因为如果您指定所需的精度,那么根据 weibull 的 cdf,您可以轻松估计尾部的多少截短。如果你有一个固定的限制,那么你可以准确(或几乎)指定细分的数量(例如,为了每个周期有几个(4-8)个点)。

于 2012-09-20T02:15:45.767 回答
2

我和 Jay 有同样的问题——不是 Weibull 分布,而是集成函数。我在对此问题的评论中找到了对Jay 问题 3的回答:

R中的发散积分在Wolfram中是可解的

R 包pracma包含几个用于数值求解积分的函数。在包中,可以找到一些用于集成某些数学函数的 R 函数。还有一个更一般的函数积分。这对我有帮助。示例代码如下。

问题 2:链接问题的第一个答案(上面)指出,R 打印出的不是 C 源文件的完整错误消息(该函数可能收敛速度太慢)。因此,我同意 Jay 的看法,即余弦的快速波动可能是一个问题。就我而言,在下面的示例中,这就是问题所在。

示例代码

# load Practical Numerical Math Functions package
library(pracma)

# define function
testfun <- function(r) cos(r*10^6)*exp(-r)

# Integrate it numerically with the basic 'integrate'.
out1 = integarte(testfun, 0, 100)
# "Error in integrate(testfun, 0, 100) : the integral is probably divergent"

# Integrate it numerically with 'integral' from 'pracma' package
# using 'Gauss-Kronrod' method and 10^-8 as relative tolerance. I
# did not try the other ones.
out2 = integral(testfun, 0, 100, method = 'Kronrod', reltol = 1e-8)

两个备注

  1. 积分函数不会像积分函数那样中断,它可能需要很长时间才能运行。我不知道(我也没有尝试过)用户是否可以限制迭代次数(?)。
  2. 即使积分函数最终没有错误,我也不确定结果有多正确对一个在零附近快速波动的函数进行数值积分似乎非常棘手,因为人们不知道在哪里计算了波动函数的确切值(正值是负值的两倍;接近局部最大值的正值和远离局部最大值的负值) . 我不是数字积分方面的专家,但我只是在我的数字讲座中了解了一些基本的固定步长积分方法。因此,积分中使用的自适应方法可能以某种方式解决了这个问题。
于 2015-02-26T23:55:20.083 回答
0

我正在尝试回答问题 1 和 3。话虽如此,我并没有提供任何原始代码。我做了一个谷歌搜索,希望这会有所帮助。祝你好运!

来源:http://cran.r-project.org/doc/contrib/Ricci-distributions-en.pdf (p.6)

#Script

library(ggplot2)

## sampling from a Weibull distribution with parameters shape=2.1 and scale=1.1
x.wei<-rweibull(n=200,shape=2.1,scale=1.1) 

#Weibull population with known paramters shape=2 e scale=1
x.teo<-rweibull(n=200,shape=2, scale=1) ## theorical quantiles from a

#Figure
qqplot(x.teo,x.wei,main="QQ-plot distr. Weibull") ## QQ-plot
abline(0,1) ## a 45-degree reference line is plotted
于 2012-09-17T02:28:23.997 回答
0

这有什么用吗?

http://www.sciencedirect.com/science/article/pii/S0378383907000452

Muraleedharana et al (2007) Modified Weibull distribution for maximum and significant wave height simulation and prediction, Coastal Engineering,第 54 卷,第 8 期,2007 年 8 月,第 630-638 页

摘自:“威布尔分布的特征函数推导出来。”

于 2012-09-17T03:08:00.237 回答