3 パラメーター ワイブル分布

データの読み込み

carsmall データ セットを読み込みます。このデータ セットには、1970 年代と 1980 年代初期に製造された自動車の測定値が格納されています。

load carsmall

この例では、変数 Weight に格納された車両の重量の測定値を使用します。

2 パラメーター ワイブル分布の当てはめ

最初に、2 パラメーター ワイブル分布を Weight に当てはめます。

pd = fitdist(Weight,'Weibull')

pd = 
  WeibullDistribution

  Weibull distribution
    A = 3321.64   [3157.65, 3494.15]
    B = 4.10083   [3.52497, 4.77076]

ヒストグラムとの適合をプロットします。

figure
histogram(Weight,8,'Normalization','pdf')
hold on
x = linspace(0,6000);
plot(x,pdf(pd,x),'LineWidth',2)
hold off

当てはめ後の分布プロットがヒストグラムと十分に一致していません。ヒストグラムから Weight < 1500 の領域に標本がないことがわかります。Weight から 1500 を減算し、ワイブル分布をもう一度当てはめます。

pd = fitdist(Weight-1500,'Weibull')

pd = 
  WeibullDistribution

  Weibull distribution
    A = 1711.75   [1543.58, 1898.23]
    B = 1.99963   [1.70954, 2.33895]

figure
histogram(Weight-1500,8,'Normalization','pdf')
hold on
plot(x,pdf(pd,x),'LineWidth',2)
hold off

当てはめ後の分布プロットが前よりも一致しています。

分布の範囲に任意の値を指定する代わりに、3 パラメーター ワイブル分布のカスタム関数を定義して範囲 (位置パラメーター $$c$$) を推定できます。

3 パラメーター ワイブル分布のカスタム pdf の定義

3 パラメーター ワイブル分布の確率密度関数を定義します。

f_def = @(x,a,b,c) (x>c).*(b/a).*(((x-c)/a).^(b-1)).*exp(-((x-c)/a).^b);

あるいは、関数 wblpdf を使用して 3 パラメーター ワイブル分布を定義できます。

f = @(x,a,b,c) wblpdf(x-c,a,b);

3 パラメーター ワイブル分布の当てはめ

関数 mle を使用して、3 つのパラメーターの MLE を求めます。カスタム分布には初期パラメーター値 (名前と値の引数 Start) も指定する必要があります。

try
    mle(Weight,'pdf',f,'Start',[1700 2 1500])
catch ME
    disp(ME)
end

  MException with properties:

    identifier: 'stats:mle:NonpositivePdfVal'
       message: 'Custom probability function returned negative or zero values.'
         cause: {}
         stack: [12x1 struct]
    Correction: []

カスタム関数が非正の値を返すため、mle からエラーが返されます。このエラーは、分布の下限を当てはめたときや確率密度がゼロになる領域がある分布を当てはめたときに発生する典型的な問題です。mle は、密度がゼロになるパラメーターの値をいくつか試した後、パラメーターの推定に失敗します。上記の関数の呼び出しでは、mle は Weight の最小値よりも大きい $$c$$ の値を試し、その一部の点で密度がゼロになるためエラーを返しています。

この問題を回避するため、関数 mle を呼び出すときに、無効な関数値をチェックするオプションをオフにして、パラメーターの境界を指定できます。

関数 mle の反復推定処理に関する既定のオプションを表示します。

statset('mlecustom')

ans = struct with fields:
          Display: 'off'
      MaxFunEvals: 400
          MaxIter: 200
           TolBnd: 1.0000e-06
           TolFun: 1.0000e-06
       TolTypeFun: []
             TolX: 1.0000e-06
         TolTypeX: []
          GradObj: 'off'
         Jacobian: []
        DerivStep: 6.0555e-06
      FunValCheck: 'on'
           Robust: []
     RobustWgtFun: []
           WgtFun: []
             Tune: []
      UseParallel: []
    UseSubstreams: []
          Streams: {}
        OutputFcn: []

関数 statset で作成したオプション構造体を使用して既定値をオーバーライドします。FunValCheck フィールドの値を 'off' と指定して、カスタム関数の値についての有効性チェックをオフにします。

opt = statset('FunValCheck','off');

3 つのパラメーターの MLE をもう一度求めます。反復処理オプション (名前と値の引数 Options) とパラメーターの境界 (名前と値の引数 LowerBound および UpperBound) を指定します。スケール パラメーターと形状パラメーターは正である必要があります。位置パラメーターは標本データの最小値より小さい必要があります。

params = mle(Weight,'pdf',f,'Start',[1700 2 1500],'Options',opt, ...
    'LowerBound',[0 0 -Inf],'UpperBound',[Inf Inf min(Weight)])

params = 1×3
103 ×

    1.3874    0.0015    1.7581

ヒストグラムとの適合をプロットします。

figure
histogram(Weight,8,'Normalization','pdf')
hold on
plot(x,f(x,params(1),params(2),params(3)),'LineWidth',2)
hold off

当てはめ後の分布プロットが十分に一致しています。

$$b<1$$ の場合の 3 パラメーター ワイブル分布の当てはめ

スケール パラメーター $$b$$ が 1 より小さい場合、ワイブル分布の pdf は $$c$$ (位置パラメーター) の下限近くで無限大に近づきます。この問題は、該当する場合は区間打ち切りデータを指定することで回避できます。

cities データ セットを読み込みます。このデータには、アメリカ合衆国の 329 の都市について、生活の質を表す 9 つの異なる指標 (気候、住宅、健康、犯罪、交通、教育、芸術、レクリエーション、経済) の評価点が含まれています。どの指標でも、評価点が高いほど優れています。

load cities

7 番目のインジケーター (芸術) の MLE を求めます。

Y = ratings(:,7);
params1 = mle(Y,'pdf',f,'Start',[median(Y) 1 0],'Options',opt)

Warning: Maximum likelihood estimation did not converge.  Iteration limit exceeded.

params1 = 1×3
103 ×

    2.7584    0.0008    0.0520

警告メッセージは、推定が収束しないことを示しています。推定オプションを変更し、MLE をもう一度求めます。最大反復回数 (MaxIter) と目的関数評価の最大回数 (MaxFunEvals) を増やします。

opt.MaxIter = 1e3;
opt.MaxFunEvals = 1e3;
params2 = mle(Y,'pdf',f,'Start',params1,'Options',opt)

Warning: Maximum likelihood estimation did not converge.  Function evaluation limit exceeded.

params2 = 1×3
103 ×

    2.7407    0.0008    0.0520

pdf が下限近くで無限大に近づくため、反復はまだ収束しません。

Y の指標は、最も近い整数に丸められた値であると仮定します。すると、Y の値を区間打ち切り観測値として扱えます。Y の観測値 y は、実際の評価点が y–0.5 と y+0.5 の間にあることを示しています。各行が Y にあるそれぞれの整数を囲む区間を表している行列を作成します。

intervalY = [Y-0.5, Y+0.5];

intervalY を使用して MLE をもう一度求めます。カスタム分布を打ち切られたデータ セットに当てはめるには、pdf および cdf の両方を関数 mle に渡す必要があります。

F = @(x,a,b,c) wblcdf(x-c,a,b);
params = mle(intervalY,'pdf',f,'cdf',F,'Start',params2,'Options',opt)

params = 1×3
103 ×

    2.7949    0.0008    0.0515

収束の問題は発生せずに、関数で MLE が求められます。この当てはめは区間への確率の当てはめに基づくため、単一の点で密度が無限大に近づく問題は発生しません。この手法は、データの区間打ち切りバージョンへの変換が適した場合にのみ使用できます。

結果をプロットします。

figure
histogram(Y,'Normalization','pdf')
hold on
x = linspace(0,max(Y));
plot(x,f(x,params(1),params(2),params(3)),'LineWidth',2)
hold off

当てはめ後の分布プロットが十分に一致しています。