カスタム分布当てはめ時の数値問題の回避

適切な初期パラメーター値の指定

打ち切りデータの値は正確に把握できないため、最尤推定法にはさらに情報が必要です。特に、確率密度関数 (pdf)、累積分布関数 (cdf)、および適切な初期パラメーター値が対数尤度の計算には必要となります。関数 evpdf および evcdf を使用して、pdf および cdf を指定できます。

最初に、打ち切られていない極値データを生成します。

rng(0,'twister');
n = 50;
mu = 5;
sigma = 2.5;
x = evrnd(mu,sigma,n,1);

次に、事前設定されたカットオフ値よりも大きいすべての値を、カットオフ値に置き換えることによって、打ち切ります。このタイプの打ち切りは、タイプ II の打ち切りとして知られています。

c = (x > 7);
x(c) = 7;

打ち切られた観測値の割合を確認します。

sum(c)/length(c)

ans = 0.1200

オリジナル データの 12% が、7 のカットオフ値で右側打ち切りされています。

データのヒストグラム、たとえば棒グラフをプロットして、打ち切られた観測値を表現します。

[uncensCnts,Edges] = histcounts(x(~c),10);
censCnts = histcounts(x(c),Edges);
bar(Edges(1:end-1)+diff(Edges)/2,[uncensCnts' censCnts'],'stacked')
legend('Fully observed data','Censored data','Location','northwest')

データには打ち切られた観測値が含まれますが、打ち切られた観測値の比率は比較的小さいものです。そのため、モーメント法がパラメーター推定のための合理的な開始点となります。打ち切られていないデータの観測平均値および標準偏差に対応する mu と sigma の初期パラメーター値を計算します。

sigma0 = sqrt(6)*std(x(~c))./pi

sigma0 = 2.3495

mu0 = mean(x(~c))-psi(1).*sigma0

mu0 = 3.5629

2 つの極値分布パラメーターの最尤推定値と、約 95% の信頼区間を求めます。打ち切りベクトル、pdf、cdf、および初期パラメーター値を指定します。sigma (スケール パラメーター) は正でなければならないため、パラメーターの下限も指定する必要があります。

[paramEsts,paramCIs] = mle(x,'Censoring',c, ...
    'pdf',@evpdf,'cdf',@evcdf, ...
    'Start',[mu0 sigma0],'LowerBound',[-Inf,0])

paramEsts = 1×2

    4.5530    3.0215

paramCIs = 2×2

    3.6455    2.2937
    5.4605    3.7494

logpdf および logsf の指定

カスタム分布の当てはめには、パラメーターの初期推定が必要です。開始点が前提としてどの程度優れているかを判断するのは困難な場合があります。最尤推定とかけ離れた開始点を指定した場合、観測値のいくつかが、その開始点に対応する極値分布の裾部から大きくずれる可能性があります。そのような場合、次の状況のいずれかが発生する可能性があります。

cdf 値も非常に小さくなってアンダーフローする場合がありますが、この状況は問題になりません。

いずれの状況も、mle が対数尤度を計算するときに問題の原因となることがあります。その理由は、いずれも対数尤度値が —Inf となり、この値を mle の最適化アルゴリズムが処理できないからです。

別の開始点ではどうなるか調べます。

start = [1 1];
try
    [paramEsts,paramCIs] = mle(x,'Censoring',c, ...
        'pdf',@evpdf,'cdf',@evcdf, ...
        'Start',start,'LowerBound',[-Inf,0])
catch ME
    disp(ME.message)
end

Custom cumulative distribution function returned values greater than or equal to 1.

この場合、2 番目の問題の状況が発生しています。初期パラメーター推定において一部の cdf 値が厳密に 1 なので、対数尤度が無限大になります。名前と値の引数 Options を使用して、FunValCheck 制御パラメーターを off に設定してみることができます。off オプションは、非有限な尤度値のチェックを無効にします。しかし、この数値問題を解決する最善の方法はその根元にあります。

極値 cdf の形式は次のとおりです。

  p = 1 - exp(-exp((x-mu)./sigma))

打ち切り観測値の対数尤度への寄与は、その生存時間関数 (SF) 値の対数、すなわち log(1-cdf) です。極値分布の場合、SF の対数は -exp((x-mu)./sigma) です。log(1-(1-exp(logSF))) を計算するのではなく、対数生存時間関数を直接使用して対数尤度を計算できる場合は、cdf での四捨五入の問題を回避できます。倍精度で 1 と区別できない cdf 値となる観測値は、非ゼロ値として容易に表現できる対数生存時間関数値をもっています。たとえば、cdf 値 (1-1e-20) は倍精度で 1 に丸められます。それは、倍精度 eps が約 2e-16 であるためです。

SFval = 1e-20;
cdfval = 1 - SFval

cdfval = 1

対応する生存時間関数値の対数は容易に表現できます。

log(SFval)

ans = -46.0517

同じ状況は対数 pdf の場合にも成り立ちます。対数尤度に対する、打ち切られていない観測値の寄与は、その pdf 値の対数です。log(exp(logpdf)) を計算する代わりに、対数 pdf を直接使用すれば、倍精度で pdf 値が 0 と区別できないアンダーフロー問題を回避できます。対数 pdf は有限の負数として容易に表現可能です。たとえば、pdf 値 1e-400 は倍精度でアンダーフローします。倍精度 realmin が約 2e-308 であるためです。

logpdfval = -921;
pdfval = exp(logpdfval)

pdfval = 0

関数 mle を使用して、名前と値の引数 logpdf および logsf を設定することにより、(pdf や cdf ではなく) 対数 pdf および対数 SF でカスタム分布を指定できます。関数 pdf および関数 cdf とは異なり、対数 pdf および対数 SF には組み込み関数がありません。したがって、これらの値を計算する無名関数を作成する必要があります。

evlogpdf = @(x,mu,sigma) ((x-mu)./sigma - exp((x-mu)./sigma)) - log(sigma);
evlogsf = @(x,mu,sigma) -exp((x-mu)./sigma);

同じ開始点を使用して、極値分布の代替対数 pdf および対数 SF を指定することで、問題が解決可能になります。

start = [1 1];
[paramEsts,paramCIs] = mle(x,'Censoring',c, ...
    'logpdf',evlogpdf,'logsf',evlogsf, ...
    'Start',start,'LowerBound',[-Inf,0])

paramEsts = 1×2

    4.5530    3.0215

paramCIs = 2×2

    3.6455    2.2937
    5.4605    3.7494

この手順で、適切でない開始点の問題を必ず解決できるわけではありません。開始点は慎重に選択してください。

最適化関数 fmincon への勾配の提供

既定の設定では、mle は関数 fminsearch を使用してカスタム分布の対数尤度を最大化するパラメーター値を計算します。fminsearch は導関数を使用しない最適化アルゴリズムを使用するので、このタイプの問題に対する適切な選択肢となります。ただし、問題によっては、対数尤度関数の導関数を使用する最適化アルゴリズムを選択すると、特に開始点が最終的な回答とかけ離れている場合、最尤推定値に収束するか、最尤推定値に収束しないかに大きな差が出てきます。導関数を指定すると、収束を促進することもできます。

名前と値の引数 OptimFun を mle に fmincon として指定すると、関数 fmincon (Optimization Toolbox が必要) を使用できます。関数 fmincon には、導関数情報を使用できる最適化アルゴリズムが含まれています。fmincon でアルゴリズムを利用するには、対数尤度だけでなくその勾配を返すために記述されている対数尤度関数を使用して、カスタム分布を指定します。対数尤度関数の勾配は、そのパラメーターに関連する偏導関数のベクトルです。

この手法では、対数尤度とその勾配の両方を計算するコードを記述するための追加準備が必要になります。helper_evnegloglike という名前の関数を別のファイルに定義します。

function [nll,ngrad] = helper_evnegloglike(params,x,cens,freq)
%HELPER_EVNEGLOGLIKE Negative log-likelihood for the extreme value
% distribution.
% This function supports only the example Avoid Numerical Issues When
% Fitting Custom Distributions (customdist2demo.mlx) and might change in 
% a future release.


if numel(params)~=2
    error(message('stats:probdists:WrongParameterLength',2));
end
mu = params(1);
sigma = params(2);
nunc = sum(1-cens);
z = (x - mu) ./ sigma;
expz = exp(z);
nll = sum(expz) - sum(z(~cens)) + nunc.*log(sigma);
if nargout > 1
    ngrad = [-sum(expz)./sigma + nunc./sigma, ...
        -sum(z.*expz)./sigma + sum(z(~cens))./sigma + nunc./sigma];
end

関数 helper_evnegloglike は、対数尤度値と勾配値の両方について負数を返します。それは、mle によって、負の対数尤度が最小化するためです。

勾配に基づく最適化アルゴリズムを使用して最尤推定値を計算するため、名前と値の引数 nloglf、OptimFun、および Options を指定します。nloglf は負の対数尤度を計算するカスタム関数を指定します。OptimFun は fmincon を最適化関数として指定します。Options は、fmincon がカスタム関数の 2 番目の出力を勾配に使用することを指定します。

start = [1 1];
[paramEsts,paramCIs] = mle(x,'Censoring',c,'nloglf',@helper_evnegloglike, ...
    'Start',start,'LowerBound',[-Inf,0], ...
    'OptimFun','fmincon','Options',statset('GradObj','on'))

paramEsts = 1×2

    4.5530    3.0215

paramCIs = 2×2

    3.6455    2.2937
    5.4605    3.7493