このprobパッケージは、ベース R 分布の特性関数を数値的に評価します。ほとんどすべてのディストリビューションには、既存の式があります。ただし、いくつかのケースでは、閉じた形式のソリューションは知られていません。適切な例: ワイブル分布 (ただし、以下を参照)。
ワイブル特性関数については、基本的に 2 つの積分を計算し、それらをまとめます。
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公開された結果をメイプルで使い、brute-force-integrate-by-the-definition-with-Rメイプルのお尻を蹴って正確な答えを計算した。つまり、私は同じ答え (数値精度まで) を数分の 1 秒とさらにわずかな価格で得ることができます。
編集：
探している正確な積分を書き留めるつもりでしたが、この特定のサイトは MathJAX をサポートしていないようですので、代わりにリンクを提供します。合理的な入力(それが何を意味するにせよ)に対するワイブル分布の特性関数を数値的に評価しようとしています。値は複素数ですが、実数部分と虚数部分に分割できます。それが私が呼んでいたものです。tRpIp
最後に 1 つのコメント: ウィキペディアにはワイブル cf の式 (無限級数) がリストされており、その式は私が上で参照した論文で証明されたものと一致しますが、その級数は に対してのみ成り立つことが証明されていshape > 1ます。ケース0 < shape < 1はまだ未解決の問題です。詳細は論文を参照してください。