mlecov

最尤推定量の漸近共分散

構文

acov = mlecov(params,data,'pdf',pdf)

acov = mlecov(params,data,'logpdf',logpdf)

acov = mlecov(params,data,'nloglf',nloglf)

acov = mlecov(___,Name,Value)

説明

acov = mlecov(params,data,'pdf',pdf) は、カスタム確率密度関数 pdf で指定された分布に対して、パラメーターの最尤推定量の漸近共分散行列の近似を返します。出力 acov は p 行 p の行列です。ここで p は params のパラメーターの数です。

mlecov は、観測された data から、最尤推定 params での対数尤度のヘッシアンへの有限差分近似を計算し、そのヘッセ行列の負の逆数を返します。

acov = mlecov(params,data,'logpdf',logpdf) は、カスタムな対数確率密度関数 logpdf によって指定された分布について acov を返します。

acov = mlecov(params,data,'nloglf',nloglf) は、カスタムな負の対数尤度関数 nloglf によって指定された分布について acov を返します。

acov = mlecov(___,Name,Value) では、前の構文におけるいずれかの入力引数の組み合わせに加えて、1 つ以上の名前と値の引数を使用してオプションを指定します。たとえば、打ち切られたデータや観測値の頻度を指定できます。

例

カスタム確率密度関数 (pdf)

標本データを読み込みます。

load carbig

ベクトル Weight には 406 台の自動車の重量が含まれます。

対数正規分布の pdf を返すカスタム関数を定義します。現在のフォルダーにファイルを lognormpdf.m として保存します。

function newpdf = lognormpdf(data,mu,sigma)
newpdf = exp((-(log(data)-mu).^2)/(2*sigma^2))./(data*sigma*sqrt(2*pi));

カスタム分布のパラメーター mu および sigma を推定します。

[phat,pci] = mle(Weight,'pdf',@lognormpdf,'Start',[4.5 0.3])

phat = 1×2

    7.9600    0.2804

pci = 2×2

    7.9327    0.2611
    7.9872    0.2997

パラメーター推定の近似共分散行列を計算します。

acov = mlecov(phat,Weight,'pdf',@lognormpdf)

acov = 2×2
10-3 ×

    0.1937   -0.0000
   -0.0000    0.0968

推定値の標準誤差を推定します。

se = sqrt(diag(acov))'

se = 1×2

    0.0139    0.0098

mu と sigma の推定値の標準誤差はそれぞれ 0.0139 と 0.0098 です。

Wald 法 (正規近似) を使用して、標準誤差 se から信頼区間 pci を再計算します。

alpha = 0.05;
probs = [alpha/2; 1-alpha/2];
pci2 = norminv(repmat(probs,1,numel(phat)),[phat; phat],[se; se])

pci2 = 2×2

    7.9327    0.2611
    7.9872    0.2997

カスタム対数確率密度関数 (pdf)

ベータ分布の対数 pdf を返すカスタム関数を定義します。現在のフォルダーにファイルを betalogpdf.m として保存します。

function logpdf = betalogpdf(x,a,b)
logpdf = (a-1)*log(x)+(b-1)*log(1-x)-betaln(a,b);

パラメーター 1.23 および 3.45 をもつベータ分布から標本データを生成し、シミュレーションを実行したデータを使用してパラメーターを推定します。

rng('default') % For reproducibility
x = betarnd(1.23,3.45,25,1);
phat = mle(x,'Distribution','beta')

phat = 1×2

    1.1213    2.7182

パラメーター推定の近似共分散行列を計算します。

acov = mlecov(phat,x,'logpdf',@betalogpdf)

acov = 2×2

    0.0810    0.1646
    0.1646    0.6074

カスタムの負の対数尤度関数

標本データを読み込みます。

load('readmissiontimes.mat')

このデータには、100 人の患者の再入院時間を示す ReadmissionTime が含まれています。このデータは、シミュレーションされたものです。

パラメーター lambda をもつポアソン分布のカスタムな負の対数尤度関数を定義します。ここで 1/lambda は分布の平均値です。カスタム関数で打ち切り情報の logical ベクトルおよびデータ頻度の整数ベクトルを使用しない場合でも、これらの値を受け入れる関数を定義する必要があります。

custnloglf = @(lambda,data,cens,freq) ...
    - length(data)*log(lambda) + sum(lambda*data,'omitnan');

カスタム分布のパラメーターを推定し、初期パラメーター値 (名前と値の引数 Start) を指定します。

phat = mle(ReadmissionTime,'nloglf',custnloglf,'Start',0.05)

phat = 
0.1462

パラメーター推定の分散を計算します。

acov = mlecov(phat,ReadmissionTime,'nloglf',custnloglf)

acov = 
2.1374e-04

標準誤差を計算します。

sqrt(acov)

ans = 
0.0146

右側打ち切りデータの指定

標本データを読み込みます。

load('readmissiontimes.mat');

このデータには、100 人の患者の再入院時間を示す ReadmissionTime が含まれています。列ベクトル Censored には各患者の打ち切り情報が含まれ、1 は右側打ち切り観測を示し、0 は正確な再入院時間が観測されることを示します。このデータは、シミュレーションされたものです。

スケール パラメーター lambda および形状パラメーター k をもつワイブル分布のカスタム対数確率密度関数 (pdf) および対数生存時間関数を定義します。データに打ち切られた観測値が含まれている場合、対数 pdf および対数生存時間関数の両方を mle および mlecov に渡す必要があります。

custlogpdf = @(data,lambda,k) ...
    log(k) - k*log(lambda) + (k-1)*log(data) - (data/lambda).^k;
custlogsf = @(data,lambda,k) - (data/lambda).^k;

打ち切れらた標本データのカスタム分布のパラメーターを推定します。カスタム分布の初期パラメーター値 (名前と値の引数 Start) を指定します。

phat = mle(ReadmissionTime,'logpdf',custlogpdf,'logsf',custlogsf, ...
    'Start',[1,0.75],'Censoring',Censored)

phat = 1×2

    9.2090    1.4223

カスタム分布のスケール パラメーターと形状パラメーターはそれぞれ 9.2090 と 1.4223 です。

パラメーター推定の近似共分散行列を計算します。

acov = mlecov(phat,ReadmissionTime, ...
    'logpdf',custlogpdf,'logsf',custlogsf,'Censoring',Censored)

acov = 2×2

    0.5653    0.0102
    0.0102    0.0163

入力引数

params — パラメーター推定

ベクトル

パラメーター推定。ベクトルとして指定します。これらのパラメーター推定は最尤推定でなければなりません。たとえば、mle で返されるパラメーター推定を指定できます。

データ型: single | double

data — 標本データおよび打ち切り情報

ベクトル | 2 列の行列

分布パラメーター params の推定に使用する標本データおよび打ち切り情報。標本データのベクトルとして、または標本データおよび打ち切り情報からなる 2 列の行列として指定します。

引数 data または名前と値の引数 Censoring のいずれかを使用して、標本データの打ち切り情報を指定できます。data が 2 列の行列の場合、mlecov は引数 Censoring の値を無視します。

data の観測値の打ち切りタイプに応じて、data をベクトルまたは 2 列の行列として指定します。

• 完全に観測されたデータ — data を標本データのベクトルとして指定します。

• 完全に観測された観測値、左側打ち切り観測値、または右側打ち切り観測値を含むデータ — 標本データのベクトルとして data を指定し、各観測値の打ち切り情報を含むベクトルとして名前と値の引数 Censoring を指定します。Censoring ベクトルには 0、-1、および 1 を含められます。それぞれ完全に観測された観測値、左側打ち切り観測値、右側打ち切り観測値を参照します。

• 区間打ち切り観測値を含むデータ — data を、標本データおよび打ち切り情報からなる 2 列の行列として指定します。data の各行は、各観測値に対する可能な生存時間または故障時間の範囲を指定します。以下の値のいずれかにできます。

  • [t,t] — t で完全に観測

  • [–Inf,t] — t で左側打ち切り

  • [t,Inf] — t で右側打ち切り

  • [t1,t2] — [t1,t2] 間で区間打ち切り (t1 < t2)

mlecov は、data 内の NaN 値を無視します。また、打ち切りベクトル (Censoring) または頻度ベクトル (Frequency) に NaN 値が含まれている場合、mlecov は data 内の対応する行を無視します。

データ型: single | double

pdf — カスタム確率密度関数

関数ハンドル | cell 配列

カスタム確率分布関数 (pdf)。関数ハンドルとして、または関数ハンドルおよび関数への追加の引数を含む cell 配列として指定します。

このカスタム関数は、標本データを含むベクトル、1 つ以上の個別の分布パラメーター、および cell 配列で渡される任意の追加引数を入力パラメーターとして受け入れます。この関数は確率密度値のベクトルを返します。

例: @newpdf

データ型: function_handle | cell

logpdf — カスタム対数確率密度関数

関数ハンドル | cell 配列

カスタム対数確率密度関数。関数ハンドルとして、または関数ハンドルおよび関数への追加の引数を含む cell 配列として指定します。

このカスタム関数は、標本データを含むベクトル、1 つ以上の個別の分布パラメーター、および cell 配列で渡される任意の追加引数を入力パラメーターとして受け入れます。この関数は対数確率値のベクトルを返します。

例: @customlogpdf

データ型: function_handle | cell

nloglf — カスタムの負の対数尤度関数

関数ハンドル | cell 配列

カスタムの負の対数尤度関数。関数ハンドルとして、または関数ハンドルおよび関数への追加の引数を含む cell 配列として指定します。

このカスタム関数は、次の入力引数を表にリストされた順番で受け入れます。

カスタム関数の入力引数	説明
params	分布パラメーター値 params のベクトル。
data	標本データ。data の値は、標本データのベクトル、または標本データおよび打ち切り情報からなる 2 列の行列です。
cens	打ち切り情報の logical ベクトル。nloglf は、名前と値の引数 Censoring を使用していない場合でも、cens を受け入れなければなりません。この場合、cens を無視するように nloglf を作成できます。
freq	データの頻度の整数ベクトル。nloglf は、名前と値の引数 Frequency を使用していない場合でも、freq を受け入れなければなりません。この場合、freq を無視するように nloglf を作成できます。
trunc	打ち切り範囲の 2 要素の数値ベクトル。nloglf は、名前と値の引数 TruncationBounds を使用している場合、trunc を受け入れなければなりません。

nloglf は、cell 配列で渡される追加引数を入力パラメーターとしてオプションで受け入れることができます。

nloglf はスカラーの負の対数尤度値を返し、オプションで負の対数尤度勾配ベクトル (名前と値の引数 Options の GradObj フィールドを参照) を返します。

例: @negloglik

データ型: function_handle | cell

名前と値の引数

オプションの引数のペアを Name1=Value1,...,NameN=ValueN として指定します。ここで、Name は引数名で、Value は対応する値です。名前と値の引数は他の引数の後に指定しなければなりませんが、ペアの順序は重要ではありません。

R2021a より前では、名前と値をそれぞれコンマを使って区切り、Name を引用符で囲みます。

例: 'Censoring',cens,'Options',opt は、mlecov に対してベクトル cens から打ち切りデータ情報を読み取って、新しいオプション構造体 opt に従って実行するように指示します。

cdf — カスタム累積分布関数

関数ハンドル | cell 配列

カスタム累積分布関数 (cdf)。関数ハンドルとして、または関数ハンドルおよび関数への追加の引数を含む cell 配列として指定します。

このカスタム関数は、標本データを含むベクトル、1 つ以上の個別の分布パラメーター、および cell 配列で渡される任意の追加引数を入力パラメーターとして受け入れます。この関数は cdf 値のベクトルを返します。

打ち切られた観測値または切り捨てられた観測値の場合、cdf および pdf の両方を定義する必要があります。完全に観測された観測値および切り捨てられていない観測値の場合、mlecov は、cdf を使用しません。data または Censoring のいずれかを使用して打ち切り情報を指定したり、TruncationBounds を使用して打ち切りの範囲を指定したりできます。

例: 'cdf',@newcdf

データ型: function_handle | cell

logsf — カスタム対数生存時間関数

関数ハンドル | cell 配列

カスタム対数生存時間関数。関数ハンドルとして、または関数ハンドルおよび関数への追加の引数を含む cell 配列として指定します。

このカスタム関数は、標本データを含むベクトル、1 つ以上の個別の分布パラメーター、および cell 配列で渡される任意の追加引数を入力パラメーターとして受け入れます。この関数は対数生存確率値のベクトルを返します。

打ち切られた観測値または切り捨てられた観測値の場合、logsf および logpdf の両方を定義する必要があります。完全に観測された観測値および切り捨てられていない観測値の場合、mlecov は、logsf を使用しません。data または Censoring のいずれかを使用して打ち切り情報を指定したり、TruncationBounds を使用して打ち切りの範囲を指定したりできます。

例: 'logsf',@logsurvival

データ型: function_handle | cell

Censoring — 打ち切りデータのインジケーター

0 のベクトル (既定値) | 0、–1、および 1 からなるベクトル

打ち切りデータのインジケーター。0、–1 および 1 からなるベクトルとして指定します。それぞれ完全に観測された観測値、左側打ち切り観測値、右側打ち切り観測値を示します。Censoring 値の各要素は、data の対応する観測の打ち切りステータスを示します。Censoring 値は、data と同じサイズでなければなりません。既定の設定は 0 のベクトルで、すべての観測値が完全に観測されていることを示します。

この引数を使用しても、区間打ち切り観測値は指定できません。標本データに区間打ち切り観測値が含まれる場合、2 列の行列を使用して data を指定します。mlecov は、data が 2 列の行列の場合に Censoring 値を無視します。

打ち切りデータの場合、pdf と cdf、logpdf と logsf、または nloglf を使用して、カスタム分布を定義する必要があります。

mlecov は、打ち切りベクトルにある任意の NaN 値を無視します。さらに、data または頻度ベクトル (Frequency) に NaN 値が含まれている場合、mlecov は打ち切りベクトル内の対応する値を無視します。

例: 'Censoring',censored。ここで、censored は打ち切り情報を含むベクトルです。

データ型: logical | single | double

TruncationBounds — 打ち切り範囲

[-Inf,Inf] (既定値) | 2 要素のベクトル

打ち切り範囲。2 要素のベクトルを指定します。

打ち切りデータの場合、pdf と cdf、logpdf と logsf、または nloglf を使用して、カスタム分布を定義する必要があります。

例: 'TruncationBounds',[0,10]

データ型: single | double

Frequency — 観測の頻度

1 から成るベクトル (既定値) | 非負の整数カウントのベクトル

観測の頻度。data と同じ行数をもつ非負の整数カウントのベクトルとして指定します。Frequency 値の j 番目の要素は、data の j 番目の行が観測された回数を示します。既定値は 1 のベクトルで、data の行あたり 1 回の観測を意味します。

mlecov は、この頻度ベクトルの NaN 値をすべて無視します。さらに、data または打ち切りベクトル (Censoring) に NaN 値が含まれている場合、mlecov は頻度ベクトル内の対応する値を無視します。

例: 'Frequency',freq。ここで、freq は観測頻度を含むベクトルです。

データ型: single | double

Options — 数値オプション

statset('mlecov') (既定値) | 構造体

有限差分ヘッセ計算の数値オプション。statset で返される構造体として指定します。

関数 mlecov は、以下の statset オプションを解釈します。

フィールド名	説明
GradObj	nloglf 入力引数で与えられた関数が 2 番目の出力として負の対数尤度の勾配ベクトルを返すことができるかどうかを示すフラグ。'on' または 'off' (既定の設定) として指定します。
DerivStep	ヘッセ計算の有限差分で使われる相対ステップ サイズ。params と同じサイズのベクトルとして指定します。 既定値は eps^(1/4) です。'GradObj' が 'on' の場合、既定値よりも小さい値が適切です。

例: 'Options',statset('GradObj','on')

データ型: struct

詳細

打ち切りのタイプ

mlecov は、左側打ち切り観測値、右側打ち切り観測値、および区間打ち切り観測値をサポートします。

• 時間 t での左側打ち切り観測 — 時間 t より前に発生したイベント。厳密なイベント時間は不明です。

• 時間 t での右側打ち切り観測 — 時間 t より後に発生したイベント。厳密なイベント時間は不明です。

• 区間 [t1,t2] 内の区間打ち切り観測 — 時間 t1 より後かつ時間 t2 より前に発生したイベント。厳密なイベント時間は不明です。

二重打ち切りデータには、左側打ち切り観測値および右側打ち切り観測値の両方が含まれます。

生存時間関数

生存時間関数とは、時間の関数としての生存の確率です。これは生存時間関数とも呼ばれます。

生存時間関数は、個体の生存時間が特定の値を超える確率を示します。累積分布関数 F(t) は、生存時間が特定の時点 t 以下である確率のため、連続分布の生存時間関数 S(t) は、累積分布関数の補数となります。S(t) = 1 – F(t)。

バージョン履歴

R2006a より前に導入

参考

mle

トピック

• 最尤推定法
• 生存時間分析とは
• カスタム分布の当てはめ
• カスタム分布当てはめ時の数値問題の回避