データセットがあり、どの分布が自分のデータに最も適しているかを調べたいと思います。

fitdistr()関数を使用して、仮定された分布を記述するために必要なパラメーターを推定しました(つまり、ワイブル、コーシー、ノーマル)。これらのパラメーターを使用して、コルモゴロフ-スミルノフ検定を実行して、サンプルデータが仮定した分布と同じ分布からのものであるかどうかを推定できます。

p値が> 0.05の場合、サンプルデータは次のようになります。同じ分布から描画されます。しかし、p値は適合の神性についての情報を提供しませんね。

では、サンプルデータのp値が正規分布とワイブル分布で> 0.05の場合、どの分布がデータに適しているかをどのように知ることができますか?

これは基本的に私が行ったことです:

> mydata [1] 37.50 46.79 48.30 46.04 43.40 39.25 38.49 49.51 40.38 36.98 40.00 [12] 38.49 37.74 47.92 44.53 44.91 44.91 40.00 41.51 47.92 36.98 43.40 [23] 42.26 41.89 38.87 43.02 39.25 40.38 42.64 36.98 44.15 44.91 43.40 [34] 49.81 38.87 40.00 52.45 53.13 47.92 52.45 44.91 29.54 27.13 35.60 [45] 45.34 43.37 54.15 42.77 42.88 44.26 27.14 39.31 24.80 16.62 30.30 [56] 36.39 28.60 28.53 35.84 31.10 34.55 52.65 48.81 43.42 52.49 38.00 [67] 38.65 34.54 37.70 38.11 43.05 29.95 32.48 24.63 35.33 41.34 # estimate shape and scale to perform KS-test for weibull distribution > fitdistr(mydata, "weibull") shape scale 6.4632971 43.2474500 ( 0.5800149) ( 0.8073102) # KS-test for weibull distribution > ks.test(mydata, "pweibull", scale=43.2474500, shape=6.4632971) One-sample Kolmogorov-Smirnov test data: mydata D = 0.0686, p-value = 0.8669 alternative hypothesis: two-sided # KS-test for normal distribution > ks.test(mydata, "pnorm", mean=mean(mydata), sd=sd(mydata)) One-sample Kolmogorov-Smirnov test data: mydata D = 0.0912, p-value = 0.5522 alternative hypothesis: two-sided 

p値はワイブル分布では0.8669、正規分布。したがって、私のデータはワイブル分布と正規分布に従っていると推測できます。しかし、どの分布関数が私のデータをよりよく説明していますか?


11dollar を参照すると、次のコードが見つかりましたが、結果の解釈方法がわかりません:

fits <- list(no = fitdistr(mydata, "normal"), we = fitdistr(mydata, "weibull")) sapply(fits, function(i) i$loglik) no we -259.6540 -257.9268 

コメント

  • データに最も適した分布を見つけたいのはなぜですか?
  • 疑似生成したいので-指定された分布に続く乱数。
  • ' KSを使用して、データセットから見つかったパラメーターを持つ分布がデータセットと一致するかどうかを確認することはできません。たとえば、このページに加えて、代替手段(および、KSテストが誤解を招く可能性のある他の方法)。
  • 別のディスカッションここにサンプルからパラメータを推定するときにKSテストを適用する方法のコードサンプルがあります。
  • I used the fitdistr() function .. .. .. 'のfitdistr関数とは何ですか?Excelからのものですか?それともCで自分で書いたものですか?

回答

まず、簡単なコメントをいくつか示します。

  • $ p $ -コルモゴロフの値-推定パラメータを使用したSmirnov-Test(KS-Test)はかなり間違っています。そのため、残念ながら、分布を適合させてから、コルモゴロフ-スミルノフ検定で推定されたパラメータを使用してサンプルをテストすることはできません。
  • サンプルが特定の条件に従うことは決してありません。したがって、たとえKS-Testの $ p $ 値が有効であり、 $ > 0.05 $ は、データがこの特定の分布に従っていることを除外できないことを意味します。別の定式化は、サンプルが特定の分布と互換性があるというものです。しかし、"私のデータは分布xyに正確に従っていますか?"は常にいいえです。
  • ここでの目標は、サンプルがどの分布に従うかを確実に判断することではありません。目標は、@ whuber(コメント内)がデータの簡潔なおおよその説明と呼ぶものです。特定のパラメトリック分布を持つことは、データのモデルとして役立ちます。

しかし、いくつかの調査を行いましょう。優れた fitdistrplus パッケージは、分布フィッティングのためのいくつかの優れた関数を提供します。関数descdistを使用して取得します可能な候補分布に関するいくつかのアイデア。

library(fitdistrplus) library(logspline) x <- c(37.50,46.79,48.30,46.04,43.40,39.25,38.49,49.51,40.38,36.98,40.00, 38.49,37.74,47.92,44.53,44.91,44.91,40.00,41.51,47.92,36.98,43.40, 42.26,41.89,38.87,43.02,39.25,40.38,42.64,36.98,44.15,44.91,43.40, 49.81,38.87,40.00,52.45,53.13,47.92,52.45,44.91,29.54,27.13,35.60, 45.34,43.37,54.15,42.77,42.88,44.26,27.14,39.31,24.80,16.62,30.30, 36.39,28.60,28.53,35.84,31.10,34.55,52.65,48.81,43.42,52.49,38.00, 38.65,34.54,37.70,38.11,43.05,29.95,32.48,24.63,35.33,41.34) 

ここで、descdistを使用しましょう:

descdist(x, discrete = FALSE) 

Descdist

サンプルの尖度と歪度の二乗は、"観測"。可能な分布には、Weibull、Lognormal、場合によってはガンマ分布が含まれるようです。

ワイブル分布と正規分布:

fit.weibull <- fitdist(x, "weibull") fit.norm <- fitdist(x, "norm") 

次に、正規分布の適合度を調べます:

plot(fit.norm) 

または不適合

ワイブル適合の場合:

plot(fit.weibull) 

ワイブル適合

どちらも見栄えは良いですが、QQプロットで判断すると、ワイブルは、特にテールで少し良く見えるかもしれません。これに対応して、ワイブルフィットのAICは、通常のフィットと比較して低くなります。

fit.weibull$aic [1] 519.8537 fit.norm$aic [1] 523.3079 

コルモゴロフ-スミルノフ検定シミュレーション

ここで説明されている@Aksakalの手順を使用して、nullの下でのKS統計をシミュレートします。

n.sims <- 5e4 stats <- replicate(n.sims, { r <- rweibull(n = length(x) , shape= fit.weibull$estimate["shape"] , scale = fit.weibull$estimate["scale"] ) estfit.weibull <- fitdist(r, "weibull") # added to account for the estimated parameters as.numeric(ks.test(r , "pweibull" , shape= estfit.weibull$estimate["shape"] , scale = estfit.weibull$estimate["scale"])$statistic ) }) 

シミュレートされたKS統計のECDFは次のようになります。

plot(ecdf(stats), las = 1, main = "KS-test statistic simulation (CDF)", col = "darkorange", lwd = 1.7) grid() 

シミュレートされたKS統計

最後に、シミュレートされたヌル分布を使用した $ p $ 値KS統計の例は次のとおりです。

fit <- logspline(stats) 1 - plogspline(ks.test(x , "pweibull" , shape= fit.weibull$estimate["shape"] , scale = fit.weibull$estimate["scale"])$statistic , fit ) [1] 0.4889511 

これにより、サンプルがワイブル分布と互換性があるというグラフィカルな結論が確認されます。

説明どおりここでは、ブートストラップを使用して、推定ワイブルPDFまたはCDFに点ごとの信頼区間を追加できます。

xs <- seq(10, 65, len=500) true.weibull <- rweibull(1e6, shape= fit.weibull$estimate["shape"] , scale = fit.weibull$estimate["scale"]) boot.pdf <- sapply(1:1000, function(i) { xi <- sample(x, size=length(x), replace=TRUE) MLE.est <- suppressWarnings(fitdist(xi, distr="weibull")) dweibull(xs, shape=MLE.est$estimate["shape"], scale = MLE.est$estimate["scale"]) } ) boot.cdf <- sapply(1:1000, function(i) { xi <- sample(x, size=length(x), replace=TRUE) MLE.est <- suppressWarnings(fitdist(xi, distr="weibull")) pweibull(xs, shape= MLE.est$estimate["shape"], scale = MLE.est$estimate["scale"]) } ) #----------------------------------------------------------------------------- # Plot PDF #----------------------------------------------------------------------------- par(bg="white", las=1, cex=1.2) plot(xs, boot.pdf[, 1], type="l", col=rgb(.6, .6, .6, .1), ylim=range(boot.pdf), xlab="x", ylab="Probability density") for(i in 2:ncol(boot.pdf)) lines(xs, boot.pdf[, i], col=rgb(.6, .6, .6, .1)) # Add pointwise confidence bands quants <- apply(boot.pdf, 1, quantile, c(0.025, 0.5, 0.975)) min.point <- apply(boot.pdf, 1, min, na.rm=TRUE) max.point <- apply(boot.pdf, 1, max, na.rm=TRUE) lines(xs, quants[1, ], col="red", lwd=1.5, lty=2) lines(xs, quants[3, ], col="red", lwd=1.5, lty=2) lines(xs, quants[2, ], col="darkred", lwd=2) 

CI_Density

#----------------------------------------------------------------------------- # Plot CDF #----------------------------------------------------------------------------- par(bg="white", las=1, cex=1.2) plot(xs, boot.cdf[, 1], type="l", col=rgb(.6, .6, .6, .1), ylim=range(boot.cdf), xlab="x", ylab="F(x)") for(i in 2:ncol(boot.cdf)) lines(xs, boot.cdf[, i], col=rgb(.6, .6, .6, .1)) # Add pointwise confidence bands quants <- apply(boot.cdf, 1, quantile, c(0.025, 0.5, 0.975)) min.point <- apply(boot.cdf, 1, min, na.rm=TRUE) max.point <- apply(boot.cdf, 1, max, na.rm=TRUE) lines(xs, quants[1, ], col="red", lwd=1.5, lty=2) lines(xs, quants[3, ], col="red", lwd=1.5, lty=2) lines(xs, quants[2, ], col="darkred", lwd=2) #lines(xs, min.point, col="purple") #lines(xs, max.point, col="purple") 

CI_CDF


GAMLSSによる自動分布フィッティング

gamlss Rのパッケージは、さまざまなディストリビューションを試し、 best " GAIC(一般化された赤池情報量基準)に準拠。主な機能はfitDistです。この関数の重要なオプションは、試行される分布のタイプです。たとえば、type = "realline"を設定すると、実数直線全体で定義された実装済みのすべての分布が試行されますが、type = "realsplus"は、実数直線で定義された分布のみが試行されます。 。もう1つの重要なオプションは、GAICのペナルティであるパラメーター $ k $ です。以下の例では、パラメータ $ k = 2 $ を設定しました。これは、" best ディストリビューションは、従来のAICに従って選択されます。 $ k $ は、 $ \ log(n)$ など、好きなように設定できます。 BIC。

library(gamlss) library(gamlss.dist) library(gamlss.add) x <- c(37.50,46.79,48.30,46.04,43.40,39.25,38.49,49.51,40.38,36.98,40.00, 38.49,37.74,47.92,44.53,44.91,44.91,40.00,41.51,47.92,36.98,43.40, 42.26,41.89,38.87,43.02,39.25,40.38,42.64,36.98,44.15,44.91,43.40, 49.81,38.87,40.00,52.45,53.13,47.92,52.45,44.91,29.54,27.13,35.60, 45.34,43.37,54.15,42.77,42.88,44.26,27.14,39.31,24.80,16.62,30.30, 36.39,28.60,28.53,35.84,31.10,34.55,52.65,48.81,43.42,52.49,38.00, 38.65,34.54,37.70,38.11,43.05,29.95,32.48,24.63,35.33,41.34) fit <- fitDist(x, k = 2, type = "realplus", trace = FALSE, try.gamlss = TRUE) summary(fit) ******************************************************************* Family: c("WEI2", "Weibull type 2") Call: gamlssML(formula = y, family = DIST[i], data = sys.parent()) Fitting method: "nlminb" Coefficient(s): Estimate Std. Error t value Pr(>|t|) eta.mu -24.3468041 2.2141197 -10.9962 < 2.22e-16 *** eta.sigma 1.8661380 0.0892799 20.9021 < 2.22e-16 *** 

AICによると、ワイブル分布(より具体的には、WEI2、その特別なパラメーター化) )データに最適です。分布WEI2の正確なパラメーター化については、279ページのこのドキュメントで詳しく説明しています。 ワームプロット(基本的にトレンド除去されたQQプロット)の残差を確認します。

WormPlot

残差は中央の水平線に近く、95%は上部の間にあると予想されます。 95%の点ごとの信頼区間として機能する下の点線の曲線。この場合、ワームプロットは私にはうまく見え、Weibull分布が適切に適合していることを示しています。

コメント

  • +1素晴らしい分析。ただし、1つの質問です。特定の主要な分布(この場合はWeibull)との互換性について肯定的な結論を出すことで、混合分布の可能性を排除できます'の存在?または、適切な混合分析を実行し、GoFをチェックしてそのオプションを除外しますか?
  • @AleksandrBlekh混合物を除外するのに十分な力を持つことは不可能です。混合物が2つのほぼ同一の分布である場合、検出できず、1つを除くすべての成分の比率が非常に小さい場合。検出することもできません。通常(分布形式を示唆する可能性のある理論がない場合)、データの簡潔な近似記述を実現するために、パラメトリック分布に適合します。混合物はそれらのどれでもありません:それらはあまりにも多くのパラメータを必要とし、目的のためにあまりにも柔軟性があります。
  • @whuber:+1あなたの優れた説明に感謝します!
  • @Lourencoカレンとフェイのグラフを見ました。青い点はサンプルを示しています。ポイントがワイブル、対数正規、ガンマ(ワイブルとガンマの間にある)の線に近いことがわかります。これらの各分布を近似した後、関数gofstatとAICを使用して適合度統計を比較しました。 "最良の"分布を決定するための最良の方法についてのコンセンサスはありません'です。グラフィカルな方法とAICが好きです。
  • @Lourenco対数正規という意味ですか?ロジスティック分布(" + "記号)は、観測されたデータからかなり離れています。対数正規分布は、私が通常見る'候補でもあります。このチュートリアルでは、投稿を短くするために、'表示しないことを選択しました。対数正規分布は、ワイブル分布と正規分布の両方と比較して、より悪い適合を示しています。 AICは537.59であり、グラフも'見栄えがよくありません。

回答

プロットは、ほとんどの場合、データがどのように見えるかをよりよく理解するための良い方法です。あなたの場合、経験累積分布関数(ecdf)を、fitdistrから取得したパラメーターを使用して理論上のcdfに対してプロットすることをお勧めします。 ()。

データに対してこれを1回行い、信頼区間も含めました。これがggplot2()を使用して取得した画像です。

ここに画像の説明を入力してください

黒い線は経験累積分布関数です色付きの線は、最尤法を使用して取得したパラメーターを使用した、さまざまな分布からの累積分布関数です。線の形式がecdfとは異なり、線がecdfからかなり離れているため、指数分布と正規分布がデータに適切に適合していないことが簡単にわかります。残念ながら、他の分布は非常に近いです。しかし、logNormal線は黒い線に最も近いと言えます。距離の尺度(MSEなど)を使用すると、仮定を検証できます。

競合する分布が2つしかない場合(たとえば、プロットに最適と思われる分布を選択する場合)、 Likelihood-Ratio-Test どの分布がより適切であるかをテストします。

コメント

  • CrossValidatedへようこそ! (a)グラフィックの作成に使用したコード、および(b)グラフィックの読み取り方法を含めるように編集できると、回答がより役立つ場合があります。
  • そこに何がプロットされていますか?それはある種の指数プロットですか?
  • しかし、どの分布がデータに最も適しているかをどのように決定しますか?グラフィックによると、' logNormalとweibullのどちらがデータに最適かを判断できませんでした。
  • 疑似乱数ジェネレータを作成する場合は、その理由を説明します。経験累積分布関数を使用しませんか?観測された分布を超える数値を描画しますか?
  • グラフを額面どおりに見ると、候補分布のどれもデータにまったく適合していないように見えます。また、ecdfの水平方向の漸近線が0.03未満であるように見えますが、これは'意味がないため、'わかりません。そもそも実際にはecdfです。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です