一个复杂函数的数值积分

The prob package numerically evaluates characteristic functions for base R distributions. For almost all distributions there are existing formulas. For a few cases, though, no closed-form solution is known. Case in point: the Weibull distribution (but see below).

prob包对基本R分布的特征函数进行数值计算。对于几乎所有的分布都有现有的公式。但是，在一些情况下，没有封闭的解决方案是已知的。恰当的例子是:威布尔分布(见下文)。

For the Weibull characteristic function I essentially compute two integrals and put them together:

对于威布尔特征函数，我本质上是计算两个积分并把它们放在一起:

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

Yes, it's clumsy, but it works surprisingly well! ...ahem, most of the time. A user reported recently that the following breaks:

是的，它很笨拙，但是它的效果却出奇的好!…嗯，大多数时候。一名用户最近报告说，以下中断:

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

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

Now, we know that the integral isn't divergent, so it must be a numerical problem. With some fiddling I could get the following to work:

现在，我们知道积分并不是发散的，所以它一定是一个数值问题。通过一些小小的改动，我可以让以下内容发挥作用:

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

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

That's OK, but it isn't quite right, plus it takes a fair bit of fiddling which doesn't scale well. I've been investigating this for a better solution. I found a recently published "closed-form" for the characteristic function with scale > 1 (see here), but it involves Wright's generalized confluent hypergeometric function which isn't implemented in R (yet). I looked into the archives for integrate alternatives, and there's a ton of stuff out there which doesn't seem very well organized.

这是可以的，但它不是很正确，而且它需要做一些不太好伸缩的小改动。我一直在寻找更好的解决办法。我发现了一个最近出版的关于> 1尺度的特征函数的“封闭形式”(见这里)，但是它涉及到Wright的广义并汇合超几何函数，这个函数在R中还没有实现。我查了一下档案库，寻找替代方案，有很多东西看起来都不太有条理。

As part of that searching it occurred to me to translate the region of integration to a finite interval via the inverse tangent, and voila! Check it out:

作为搜索的一部分，我想到了通过切反变换把积分区域转换成一个有限的区间，瞧!检查一下:

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

Questions:

Can you do better than this?
你能做得比这更好吗?
Is there something about numerical integration routines that people who are expert about such things could shed some light on what's happening here? I have a sneaking suspicion that for large t the cosine fluctuates rapidly which causes problems...?
有什么关于数值积分例程的人对于这类事情的专家能够解释这里发生的事情吗?我偷偷地怀疑，对于大的t，余弦波动很快，这会导致问题……
Are there existing R routines/packages which are better suited for this type of problem, and could somebody point me to a well-placed position (on the mountain) to start the climb?
是否存在更适合这种类型的问题的R例程/包，是否有人能给我指出一个合适的位置(在山上)开始攀登?

Comments:

Yes, it is bad practice to use t as a function argument.
是的，用t作为函数参数是不好的做法。
I calculated the exact answer for shape > 1 using the published result with Maple, and the brute-force-integrate-by-the-definition-with-R kicked Maple's ass. That is, I get the same answer (up to numerical precision) in a small fraction of a second and an even smaller fraction of the price.
我计算出了形状> 1的确切答案，使用的是与Maple公布的结果，以及基于brute-force-integrate-by- definition- r的结合。

Edit:

I was going to write down the exact integrals I'm looking for but it seems this particular site doesn't support MathJAX so I'll give links instead. I'm looking to numerically evaluate the characteristic function of the Weibull distribution for reasonable inputs t (whatever that means). The value is a complex number but we can split it into its real and imaginary parts and that's what I was calling Rp and Ip above.

我要写下我要找的精确积分，但是这个特定的网站不支持MathJAX，所以我会给出链接。我想用数值方法计算威布尔分布的特征函数为合理的输入t。这个值是一个复数但是我们可以把它分成实部和虚部这就是我上面所说的Rp和Ip。

One final comment: Wikipedia has a formula listed (an infinite series) for the Weibull c.f. and that formula matches the one proved in the paper I referenced above, however, that series has only been proved to hold for shape > 1. The case 0 < shape < 1 is still an open problem; see the paper for details.

最后一个评论是:*为威布尔c.f.列出了一个公式(一个无穷级数)，这个公式与我在上面引用的论文中证明的公式相匹配，然而，这个级数仅被证明适用于> 1型。情形0 < shape < 1仍然是一个开放的问题;详见论文。

4 个解决方案

#1

You may be interested to look at this paper, which discuss different integration methods for highly oscillating integrals -- that's what you are essentially trying to compute: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.8.6944

您可能会对这篇文章感兴趣，它讨论了用于高振荡积分的不同集成方法——这就是您试图计算的内容:http://citeseerx.ist.psuedu./viewdoc/summary? doi=10.1.1.8.6944

Also, another possible advice, is that instead of infinite limit you may want to specify a smaller one, because if you specify the precision that you want, then based on the cdf of the weibull you can easily estimate how much of the tail you can truncate. And if you have a fixed limit, then you can specify exactly (or almost) the number of subdivisions (e.g. in order to have a few(4-8) points per period).

另外，另一个可能的建议是，您可能想要指定一个更小的限制，而不是无限的限制，因为如果您指定您想要的精度，那么基于weibull的cdf，您可以轻松地估计您可以截断多少尾。如果你有一个固定的极限，那么你可以精确地(或几乎)指定细分的数量(例如为了每个周期有几个(4-8)点)。

#2

I had the same problem than Jay - not with the Weibull distribution but with the integrate function. I found my answer to Jay's question 3 in a comment to this question:

我和Jay有同样的问题-不是威布尔分布而是积分函数。我在对这个问题的评论中找到了杰伊问题3的答案:

Divergent Integral in R is solvable in Wolfram

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

The R package pracma contains several functions for solving integrals numerically. In the package, one finds some R functions for integrating certain mathematical functions. And there is a more general function integral. That helped in my case. Example code is given below.

R包实际包含了几个函数，用于数值求解积分。在这个包中，我们找到了一些R函数用来积分某些数学函数。还有一个更一般的函数积分。这对我有帮助。示例代码如下所示。

To questions 2: The first answer to the linked question (above) states that not the complete error message of the C source file is printed out by R (The function may just converge too slowly). Therefore, I would agree with Jay that the fast fluctuation of the cosine may be a problem. In my case and in the example below it was the problem.

对于问题2:链接问题(上面)的第一个答案表明C源文件的完整错误消息不是由R打印出来的(函数的收敛速度可能太慢)。因此，我同意Jay的观点，即余弦的快速波动可能是一个问题。在我的例子中，在下面的例子中是问题。

Example Code

示例代码

# 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)

Two remarks

两个评论

The integral function does not break as the integrate function does but it may take quite a long time to run. I do not know (and I did not try) whether the user can limit the number of iterations (?).
积分函数不像积分函数那样会断裂，但它可能需要很长时间才能运行。我不知道(我也没有尝试)用户是否可以限制迭代次数(?)。
Even if the integral function finalises without errors I am not sure how correct the result is. Numerically integrating a function which is fast fluctuating around zero seems to be quite tricky since one does not know where exactly values on the fluctuating function are calculated (twice as much positive than negative values; positive values close to local maxima and negative values far off). I am not on expert on numeric integration but I just got to know some basic fixed-step integration methods in my numerics lectures. So maybe the adaptive methods used in integral deal with this problem in some way.
即使积分函数没有误差，我也不知道结果是怎样的。对一个在0附近快速波动的函数进行数值积分似乎很棘手，因为人们不知道波动函数的确切值是在哪里计算的(正值是负值的两倍);正的值接近局部极值，负的值远离局部极值)。我不是数值积分的专家，但我只是在我的数字课上了解了一些基本的固定步积分方法。也许积分中的自适应方法在某种程度上解决了这个问题。

#3

I'm attempting to answer questions 1 & 3. That being said I am not contributing any original code. I did a google search and hopefully this is helpful. Good luck!

我正试图回答问题1和3。也就是说我没有贡献任何原始代码。我做了谷歌搜索，希望这对你有帮助。好运！

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

来源:http://cran.r-project.org/doc/contrib/Ricci-distributions-en.pdf(六年级)

#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

#4

Is this of any use?

这有什么用吗?

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, Volume 54, Issue 8, August 2007, Pages 630–638

Muraleedharana et al(2007)修改了Weibull分布，最大和显著的波浪高度模拟和预测，海岸工程，第54卷，第8期，2007年8月，第630-638页。

From the abstract: "The characteristic function of the Weibull distribution is derived."

文摘:导出了威布尔分布的特征函数。

#1