主要内容

カスタム分布の当てはめ

この例では,関数大中型企业を使用してカスタム分布を一変量デ,タに当てはめる方法を示します。

関数大中型企业を使用して,最尤パラメーター推定を計算したり,その精度を組み込み分布およびカスタム分布に対して推定したりできます。カスタム分布を当てはめるには,ファイル内に,または無名関数を使用して,カスタム分布に対する関数を定義する必要があります。最も簡単なケースでは,コードを記述して,当てはめる分布の確率密度関数(pdf)またはpdfの対数を計算してから,大中型企业を呼び出して分布を当てはめることができます。この例では,pdfまたはpdfの対数を使用して,以下の場合にいて説明します。

  • 打切りデタへの分布の当てはめ

  • 2の分布から成る混合の当てはめ

  • 重み付け分布の当てはめ

  • パラメタ変換を使用した,サズの小さい標本に対するパラメタ推定の正確な信頼区間の計算

カスタム関数を定義するのではなく,打切られたデタに対する大中型企业の名前と値の引数TruncationBoundsを使用できることに注意してください。また,2つの正規分布の混合には、関数fitgmdistを使用できます。この例では,これらの場合に関数大中型企业およびカスタム関数を使用します。

0を打切ったポアソン分布の当てはめ

ほとんどの場合,カウントデ,タはポアソン分布を使用してモデル化されます。関数poissfitまたはfitdistを使用して,ポアソン分布を当てはめることができます。しかし0であるカウントはデータに記録されず,0の欠損のためにポアソン分布を当てはめるのが難しい場合があります。その場合,関数大中型企业およびカスタム分布関数を使用して,0を打切ったデタにポアソン分布を当てはめます。

最初に,ランダムなポアソンデ,タを生成します。

rng (18,“旋风”%用于再现性λ = 1.75;N = 75;X1 = poissrnd(lambda,n,1);

次に,打切りをシミュレトするデタから0をすべて削除します。

X1 = X1 (X1 > 0);

打切りの後,x1の標本数をチェックします。

长度(x1)
Ans = 65

シミュレ,ションデ,タのヒストグラムをプロットします。

直方图(x1, 0:1:马克斯(x1) + 1)

图中包含一个轴对象。坐标轴对象包含一个直方图类型的对象。

デ,タは,0がないこと以外はポアソン分布のように見えます。正の整数ではポアソン分布と一致し,0では確率のないカスタム分布を使用できます。カスタム分布を使用すると,0の欠損を考慮しながら,ポアソン分布のパラメ,タ,λを推定できます。

0を打切ったポアソン分布は,その確率質量関数(pmf)で定義する必要があります。ポアソン分布の平均のパラメタλの値を指定して,x1の各点の確率を計算する無名関数を作成します。0を打ち切ったポアソン分布の及は,和が1になるように正規化されたポアソン及です。0を打切ると,正規化は1-Probability (x1 < 0)です。

Pf_truncpoiss = @(x1,lambda) poisspdf(x1,lambda)./(1-poisscdf(0,lambda));

簡略化するために,この関数に与えられるすべてのx1値は正の整数で,チェックを行わないとします。エラーチェックの場合,または複数行のコードを要するさらに複雑な分布の場合は,関数を別のファイルで定義する必要があります。

パラメタλに対し,適切なおおよその初期推定値を見けます。ここでは,標本平均を使用します。

起始=均值(x1)
Start = 2.2154

大中型企业にデータ,カスタム及関数,初期パラメーター値,およびパラメーターの下限を提供します。ポアソン分布の平均パラメ,タ,は正でなければならないため,λの下限も指定する必要があります。関数大中型企业は,λの最尤推定と,オプションとしてパラメ,タ,の近似した95%信頼区間を出力します。

[lambdaHat,lambdaCI] = mle(x1,“pdf”pf_truncpoiss,“开始”开始,...下界的, 0)
lambdaHat = 1.8760
lambdaCI =2×11.4990 - 2.2530

パラメ,タ,の推定値は,標本平均より小さくなっています。最尤推定ではデタに存在しない0を考慮します。

代わりに,名前と値の引数TruncationBoundsを使用して,打切りの範囲を指定できます。

[lambdaHat2,lambdaCI2] = mle(x1,“分布”,“泊松”,...“TruncationBounds”[0正])
lambdaHat2 = 1.8760
lambdaCI2 =2×11.4990 - 2.2530

また,mlecovが返す大きい標本の分散の近似を使用して,λの標準誤差推定も計算できます。

avar = mlecov(lambdaHat,x1,“pdf”, pf_truncpoiss);Stderr = sqrt(avar)
Stderr = 0.1923

視覚的に当てはめをチェックするために,生データの正規化ヒストグラムに対して当てはめられた及をプロットします。

直方图(x1,“归一化”,“pdf”) xgrid = min(x1):max(x1);pmfgrid = pf_truncpoiss(xgrid,lambdaHat);持有情节(xgrid pmfgrid,“- - -”)包含(x1的) ylabel (“概率”)传说(样本数据的,的安装及,“位置”,“最佳”)举行

图中包含一个轴对象。axis对象包含2个直方图类型的对象,line。这些对象表示样本数据,拟合pmf。

上限を打切った正規分布の当てはめ

連続デタは打切られている場合があります。たとえば,デ,タ収集の制限のため,所定の固定値より大きい観測値が記録されない場合があります。

その場合,打切り処理を行った正規分布のデタのシミュレションを行います。最初に,無作為な正規デ,タを生成します。

N = 500;Mu = 1;σ = 3;rng (“默认”%用于再现性X2 = normrnd(mu,sigma,n,1);

次に,打切り点xTruncを超える観測値を削除します。xTruncは推定する必要がない既知の値であると仮定します。

xTrunc = 4;x2 = x2(x2 < xTrunc);

打切りの後,x2の標本数をチェックします。

长度(x2)
Ans = 430

シミュレ,ションデ,タのヒストグラムを作成します。

直方图(x2)

图中包含一个轴对象。坐标轴对象包含一个直方图类型的对象。

x2 < xTruncでは正規分布と一致し,xTruncを超えると0確率となるカスタム分布で,シミュレションデタを当てはめます。カスタム分布を使用すると,裾の欠損を考慮しながら,正規分布のパラメ,タ,μおよびσを推定できます。

打切り処理を行った正規分布を,そのPDFで定義します。パラメタμおよびσの値を指定して,xの各点の確率密度値を計算する無名関数を作成します。打切り点は固定かつ既知であり、打ち切られた正規分布の pdf は、打ち切った後に積分が 1 となるよう正規化された pdf です。正規化はxTruncで評価されたCDFです。簡略化するために,すべてのx2値はxTruncより小さく,チェックを行わないとします。

Pdf_truncnorm = @(x2,mu,sigma)...normpdf (x2,μ、σ)。/ normcdf (xTrunc、μ、σ);

打切り点xTruncは推定する必要がないため,カスタムPDF関数の入力分布パラメ,タ,に含まれません。また,xTruncはデ,タベクトル入力引数の一部ではありません。無名関数はワ,クスペ,ス内の変数にアクセスできるので,xTruncを無名関数に追加の引数として渡す必要がありません。

パラメ,タ,推定としておおよその初期推定値を指定します。この場合,打切りが極端でないため,標本平均と標準偏差を使用します。

Start = [mean(x2),std(x2)]
开始=1×20.1585 - 2.4125

大中型企业にデータ,カスタムpdf関数,初期パラメーター値,およびパラメーターの下限を提供します。σは正でなければならないため,パラメ,タ,の下限も指定する必要があります。大中型企业は単一ベクトルとしてμおよびσの最尤推定,およびこの2のパラメタにいて約95%の信頼区間の行列を返します。

[参数,paramCIs] = mle(x2,“pdf”pdf_truncnorm,“开始”开始,...下界的(从0))
paramEsts =1×21.1298 - 3.0884
paramCIs =2×20.5713 2.7160 1.6882 3.4607

μおよびσの推定値が標本平均および標準偏差より大きくなっています。モデルの当てはめは分布の欠損した上裾を考慮しています。

代わりに,名前と値の引数TruncationBoundsを使用して,打切りの範囲を指定できます。

[paramest2,paramCIs2] = mle(x2,“分布”,“正常”,...“TruncationBounds”(负无穷xTrunc))
paramEsts2 =1×21.1297 - 3.0884
paramCIs2 =2×20.5713 2.7160 1.6882 3.4607

mlecovを使用して,パラメ,タ,推定の近似共分散行列を計算できます。通常,近似は標本が大きい場合に良好に機能するので,標準誤差は対角要素の平方根で概算できます。

acov = mlecov(参数sts,x2,“pdf”pdf_truncnorm)
acov =2×20.0812 0.0402 0.0402 0.0361
Stderr = sqrt(diag(acov))
stderr =2×10.2849 - 0.1900

視覚的に当てはめをチェックするために,生データの正規化ヒストグラムに対して当てはめたpdfをプロットします。

直方图(x2,“归一化”,“pdf”) xgrid = linspace(min(x2),max(x2));pdfgrid = pdf_truncnorm(xgrid, parametersts (1), parametersts (2));持有情节(xgrid pdfgrid,“- - -”)包含(“x2”) ylabel (的概率密度)传说(样本数据的,“安装pdf”,“位置”,“最佳”)举行

图中包含一个轴对象。axis对象包含2个直方图类型的对象,line。这些对象代表样本数据,拟合pdf。

2の正規分布から成る混合の当てはめ

データセットには二峰性または多峰性を示すものがありますが,そのようなデータに標準分布を当てはめるのは多くの場合不適切です。ただし,単純な単峰性分布を混合することで,このようなデ,タをモデル化できることがよくあります。

。次の構造的定義をもシミュレションデタを考えます。

  • まず,偏ったコ。

  • コesc escンの表が出た場合,平均 μ 1 および標準偏差 σ 1 の正規分布から値を無作為に選択します。

  • コesc escンの裏が出た場合,平均 μ 2 および標準偏差 σ 2 の正規分布から値を無作為に選択します。

当てはめているものと同じモデルを使用するのではなく,スチュ,デントの“t”分布の混合からデ,タセットを生成します。異なる分布を使用することで,モンテカルロシミュレーションで使用される手法と同様に,当てはめられているモデルの仮定からの逸脱に対して,当てはめ手法がどの程度ロバストであるかについて検定できます。

rng (10)%用于再现性X3 = [trnd(20,1,50) trnd(4,1,100)+3];直方图(x3)

图中包含一个轴对象。坐标轴对象包含一个直方图类型的对象。

確率密度を計算する無名関数を作成することにより,当てはめるモデルを定義します。2つの正規分布を混合した pdf は 2 つの正規成分の pdf の重み付き和であり、混合の確率で重み付けられます。この無名関数には、pdf を評価するデータのベクトル、および 5 つの分布パラメーターという 6 つの入力引数があります。各成分にはその平均および標準偏差のパラメーターがあります。

pdf_normmix = @(x3,p,mu1,mu2,sigma1,sigma2)...P *normpdf(x3,mu1,sigma1) + (1-p)*normpdf(x3,mu2,sigma2);

また,パラメ,タ,の初期推定値も必要です。モデルパラメ、タ、の数が増加するに、れて、開始点の定義がより重要になります。ここでは,デ,タの2の四分位点を中心とした,等しい標準偏差をも,正規分布の等量混合(p= 0.5)から開始します。標準偏差の開始値は,各成分の平均および分散として,混合の分散式から生成されます。

pStart = .5;muStart =分位数(x3,[.]25。)
muStart =1×20.3351 - 3.3046
sigmaStart =√(var(x3) - .25*diff(muStart).^2)
sigmaStart = 1.1602
start = [pStart muStart sigmaStart sigmaStart];

混合の確率に対する0と1の境界,および標準偏差の0の下限を指定します。範囲ベクトルの残りの要素は+正またはに設定し,制限がないことを示します。

lb = [0 -Inf -Inf 0 0];ub = [1 Inf Inf Inf Inf];参数= mle(x3,“pdf”pdf_normmixture,“开始”开始,...下界的磅,“UpperBound”乌兰巴托)
警告:最大似然估计不收敛。超出迭代限制。
paramEsts =1×50.3273 -0.2263 2.9914 0.9067 1.2059

警告メッセ,ジは,この関数が既定の反復設定では収束しないことを示しています。既定のオプションを表示します。

statset (“mlecustom”
ans =带字段的结构:显示:'off' MaxFunEvals: 400 MaxIter: 200 TolBnd: 1.0000 -06 TolFun: 1.0000 -06 TolTypeFun: [] TolX: 1.0000 -06 TolTypeX: [] GradObj: 'off' Jacobian:[]衍生步骤:6.0555e-06 FunValCheck: 'on' Robust: [] RobustWgtFun: [] WgtFun: [] Tune: [] UseParallel: [] UseSubstreams: [] Streams: {} OutputFcn: []

カスタム分布に対する既定の最大反復回数は200です。関数statsetで作成したオプション構造体を使用し、この既定値をオ、バ、ラ、ドして反復回数を増やします。また,最大関数評価を増やします。

选项= statset(“麦克斯特”, 300,“MaxFunEvals”, 600);参数= mle(x3,“pdf”pdf_normmixture,“开始”开始,...下界的磅,“UpperBound”乌兰巴托,“选项”选项)
paramEsts =1×50.3273 -0.2263 2.9914 0.9067 1.2059

収束の最終的な反復は,結果の最後の数桁にしか関係ありません。ただし,収束に達していることを必ず確認することをお勧めします。

視覚的に当てはめをチェックするために,生データの確率のヒストグラムに対して当てはめた密度をプロットします。

直方图(x3,“归一化”,“pdf”)举行Xgrid = linspace(1.1*min(x3),1.1*max(x3),200);Pdfgrid = pdf_normmix (xgrid,...paramEsts paramEsts (1) (2), paramEsts (3), paramEsts (4), paramEsts (5));情节(xgrid pdfgrid,“- - -”)举行包含(“x3”) ylabel (的概率密度)传说(样本数据的,“安装pdf”,“位置”,“最佳”

图中包含一个轴对象。axis对象包含2个直方图类型的对象,line。这些对象代表样本数据,拟合pdf。

代わりに,正規分布の混合には関数fitgmdistを使用できます。初期推定および反復アルゴリズムの設定により,推定値が異なる場合があります。

Mdl = fitgmdist(x3',2)
Mdl =一维2组分的高斯混合分布组件1:混合比例:0.329180均值:-0.2200组件2:混合比例:0.670820均值:2.9975
Mdl。σ
Ans = Ans (:,:,1) = 0.8274 Ans (:,:,2) = 1.4437

精度が異なるデ,タへの重み付き正規分布の当てはめ

10個のデータポイントがあり,それぞれが1 ~ 8個の観測値のいずれかの平均であるとします。元の観測値は利用できませんが,各デ,タポ,ントの観測値の個数は既知です。各ポesc escントの精度は,対応する観測値の個数によって異なります。生デ,タの平均値と標準偏差を推定する必要があります。

X4 = [0.25 -1.24 1.38 1.39 -1.43 2.79 3.52 0.92 1.44 1.26]';M = [8 2 1 3 8 4 2 5 2 4]';

各デタポントの分散は対応する観測値の個数に反比例するため,1 / mを使用して最尤近似の各デタポントの分散を重み付けします。

W = 1./m;

このモデルでは,分布はそのPDFで定義できます。ただし,正規PDFは次の形式のため,PDFの対数を使用する方がより適切です。

C .* exp(-0.5 .* z.²),

また,大中型企业はPDFの対数を取り,対数尤度を計算します。したがって,代わりにPDFの対数を直接計算する関数を作成します。

関数PDFの対数は,正規分布のパラメ,タ,μおよびσの値を指定して,xの各点に対する確率密度の対数を計算する必要があります。また,別の分散の重みを考慮する必要もあります。helper_logpdf_wn1という名前の関数を別のファ@ @ルhelper_logpdf_wn1.mに定義します。

函数Logy = helper_logpdf_wn1(x,m,mu,sigma)权重正态分布pdf的对数此函数仅支持示例适配自定义分发金宝app% (customdist1demo.mlx),并可能在将来的版本中更改。V = sigma。^2 ./ m;Logy = -(x-mu)^2 ./ (2.*v) - .5.*log(2.*pi.*v);结束

パラメ,タ,推定のおおよその初期推定値を指定します。この場合,重み付けされていない標本平均と標準偏差を使用します。

Start = [mean(x4),std(x4)]
开始=1×21.0280 - 1.5490

σは正でなければならないため,パラメ,タ,の下限を指定する必要があります。

[paramest1,paramCIs1] = mle(x4,“logpdf”,...@ (x,μ、σ)helper_logpdf_wn1 (x, m,μ、σ),...“开始”开始,下界的(从0))
paramEsts1 =1×20.6244 - 2.8823
paramCIs1 =2×2-0.2802 1.6191 1.5290 4.1456

μの推定値は標本平均の推定値の3分の2未満です。推定値に影響を与えるのは,最も信頼性の高いデ,タ点,まり最も多い生の観測値に基づくデ,タ点です。このデータセットでは,これらの点が重み付けされていない標本平均から推定値を引き下げる傾向にあります。

パラメ,タ,変換を使用した正規分布の当てはめ

関数大中型企业は,厳密法が利用できない場合,推定器の分布について,大きい標本の正規近似を使用してパラメーターの信頼区間を計算します。標本サズが小さい場合は,1以上のパラメタを変換することにより正規近似を改良できます。この場合,正規分布のスケ,ルパラメ,タ,を対数に変換します。

最初に,変換されたパラメ,タ,をσに使用するhelper_logpdf_wn2という名前の新しい対数関数PDFを作成します。

函数Logy = helper_logpdf_wn2(x,m,mu,logsigma)的权值正态分布的对数% log(sigma)参数化此函数仅支持示例适配自定义分发金宝app% (customdist1demo.mlx),并可能在将来的版本中更改。V = exp(logsigma)^2 ./ m;Logy = -(x-mu)^2 ./ (2.*v) - .5.*log(2.*pi.*v);结束

σの新しいパラメ,タ,化に変換された同じ開始点,まり,標本標準偏差の対数を使用します。

Start = [mean(x4),log(std(x4))]
开始=1×21.0280 - 0.4376

σは任意の正の値であるため,日志(σ)は非有界であり,下限や上限を指定する必要はありません。

[paramest2,paramCIs2] = mle(x4,“logpdf”,...@ (x,μ、σ)helper_logpdf_wn2 (x, m,μ、σ),...“开始”,开始)
paramEsts2 =1×20.6244 - 1.0586
paramCIs2 =2×2-0.2802 0.6203 1.5290 1.4969

パラメ,タ,化では日志(σ)が使用されるため,変換してオリジナルのスケ,ルに戻し,σの推定値と信頼区間を取得する必要があります。

= exp(paramest2 (2))
sigmaHat = 2.8823
sigmaCI = exp(paramCIs2(:))
sigmaCI =2×11.8596 - 4.4677

最尤推定がパラメ,タ,化に対し不変であるため,μおよびσ両方の推定値は最初の当てはめと同じです。σに対する信頼区間はparamCIs1 (: 2)とは少し異なります。

参考

|

関連するトピック