2

データが特定の分布関数に適合しているかどうかを視覚的に評価したいと考えています。これを行うために、R を使用して分位 - 分位 (QQ) プロットを生成しています。分布関数は非常に特殊で、確率分布の標準リストには含まれていないため、独自の R 関数を記述して記述しました。以下のコードでは「DistFunc」と呼ばれ、2 つのガンマ関数の比率で構成されます。

簡単に言うと、私のコードで行っていることは、2 つの列を含むファイル 'DistributionEstimate.txt' からデータを読み取ることです。列 1 は x の値で、列 2 は y の値です。変数 'a' と 'b' は、この分布関数のデータへの最小二乗フィットを使用して別のプログラムで以前に決定した最適なパラメーターです。次に、DistFunc を定義し、qqmath 関数を使用して QQ プロットをプロットしようとします。

この時点で問題が発生します。R は、DistFunc が 'gammafn' の範囲外の値を返し、何もプロットできないという多くの警告を表示し続けます。関数には原点に近い極が含まれていることがわかっているので、これで十分です。コードでわかるように、DistFunc を正規化して確率分布に変換しようとしています (これは、qqmath を使用するために必要なものだと思いますか?)。しかし、これは役に立ちません。

この問題を克服する方法を知っている人はいますか?たとえば、正規化を必要としない別のプロット関数を使用したり、結果に深刻な影響を与えずに疑似確率分布に変換したりできますか?

役立つ情報を提供していただければ幸いです。

install.packages('lattice')
library(lattice)
x<-read.table("C:/DistributionEstimate.txt", colClasses = c(rep("NULL",1),rep("numeric",1)), header = FALSE)
y<-read.table("C:/DistributionEstimate.txt", colClasses = c(rep("numeric",1),rep("NULL",1)), header = FALSE)
x<-sapply(x, as.numeric)
y<-sapply(y, as.numeric)
a<-16359727025.407821410;
b<-198838619.13262583836;
DistFunc <- function(k,ampl=a,stretch=b) {
    fdist<-ampl*gamma(k*stretch-1/2)/gamma(k*stretch+1)
    fnorm<-fdist/sum(fdist)
}
qqmath(DistFunc(x), y, col="blue", envelope=.95, xlab="Quantiles of the best-fit model", ylab="Quantiles of the data")
abline(0,1, col="red", lwd=2)
grid()
4

1 に答える 1

3

QQ プロットの背後にある考え方は、特定の分布から生じると考えられる観測値を、同じサイズのサンプルでその分布から見られると予想される値と比較することです。

したがって、最初の問題は、 と の両方xのy値を持っていることです。QQ プロットは単変量プロットです。分布に対して 1 つの値のセットを照合しています。ペアをプロットするための 2 番目の次元(x,y)は、分布関数によって計算されます。

qqmath期待される分布関数は密度関数ではありません。分位数を分布からの値に変換する関数が必要です。これは、やq*のように、R で機能する分布関数のファミリと同じです。関数は 0 ～ 1 の範囲の数値を受け入れ、それを の分布の定義域または の値に変換する必要があります。プロット中に、分位数のリストをこの関数に渡し、期待値のリストを取得します。次に、(ソートされた) 観測値に対して期待値のリストをプロットします。qnromqexp(-Inf,Inf)qnorm(0, Inf)qexpqqmath

qexp例として、関数を「カスタム」分位関数として使用します。観察する

myDist<-function(x) {
    qexp(x, 5)
}

set.seed(15)
x <- rexp(100, 5)
qqmath(~x, distribution=myDist, main="qqmath")

そして、これはまったく同じです

exp.x <- myDist(ppoints(length(x)))
xyplot(sort(x)~exp.x, main="xyplot")

qqmath vc xyplot

あなたが抱えている問題の1つはDistFunc、分位関数よりも密度のように見えることだと思います. 密度関数から確率に移行するには、積分する必要があります。q-likeHere'aは、任意の密度関数の関数を作成しようとするヘルパー関数です。

getq <- function(density, from, to, steps=1000) {
    x <- seq(from=from, to=to, length.out=steps) 
    y <- mapply(function(a,b)integrate(density,a,b)$value, x[-steps], x[-1])
    approxfun(c(0,cumsum(y)),x)
}

最初のパラメーターは、1 パラメーターの密度関数です。これは、統合時に使用されます。次に、fromおよびtoパラメータは、値がゼロ以外の確率を持つ場所を指定します。次にsteps、積分を実行するポイントの数です。次にapproxfun、実際に計算したポイント数と最終q関数によって要求されたポイントの間で補間を行います。これが標準密度でどのように機能するかを見てみましょう。ここでも、指数、レート 5、密度を使用します。

myq <- getq(function(x) dexp(x,5), 0, 4)

dexp密度が 1 つのパラメーターしかとらないように、レート パラメーターでをラップする匿名関数を作成することに注意してください。ここでは 0 から 4 に進みます。その時点までに確率はほぼ 1.0 になっているからです。これで、この関数を標準のように使用できますqexp

> qexp(.5,5)
[1] 0.1386294
> myq(.5)
[1] 0.1386388

.5 に対して非常によく似た回答が得られることがわかります。それで、それはうまくいっているようです。したがって、これは、分位関数が適切な閉じた形式を持たない場合に、密度関数を分位関数に変換する簡単な方法の 1 つです。

そして最後に私が目にする問題は、あなたaのb価値観が非常に大きいということです。関数内でそれらを使用すると、gammaすぐに R が処理できない数値になります。今、あなたはgamma互いに分割しているので、それらがいくらか相殺されることを願っていますが、通常、標準バージョンを使用するとオーバーフローに遭遇します. したがって、大きな値を計算するには、対数スケールで計算し、exp()すべてが完了したら自然なスケールに戻すのがコツです。したがって、関数を次のように変更できます

DistFunc <- function(k,ampl=a,stretch=b) {
    fdist <- exp(log(ampl) + lgamma(k*stretch-1/2) - lgamma(k*stretch+1))
    fnorm <- fdist/sum(fdist)
}

lgammaは対数スケールのガンマ関数であることに注意してください。しかし、あなたのaandb値を使用しても、ほとんどの場合、それだけでは十分ではないようです。パラメータを指定して、その関数からどのように使用可能な数値を取得できるかわかりません。また、あなたの分布の範囲がどのように考えられているかもわかりません。適切な密度関数のように 1 に統合する方法が見つかりませんでした。

于 2014-05-26T01:33:26.557 に答える