Main Content

ksdensity

一変量および二変量データのカーネル平滑化関数推定値

説明

[f,xi] = ksdensity(x) は、ベクトルまたは 2 列の行列 x に格納されている標本データに対する確率密度推定値 f を返します。推定値は正規カーネル関数に基づいており、x 内のデータ範囲全体で等間隔となる点 xi において評価されます。ksdensity は、一変量データの場合は 100 個の、二変量データの場合は 900 個の点で密度を推定します。

ksdensity は、連続分布標本で最も適切に機能します。

[f,xi] = ksdensity(x,pts) は、f を評価する点 (pts) を指定します。この場合、xipts には同じ値が格納されます。

[f,xi] = ksdensity(___,Name,Value) は、前の構文の入力引数のいずれかに加えて、1 つ以上の名前と値のペアの引数によって指定される追加オプションを使用します。たとえば、確率密度関数、累積確率関数、生存時間関数のような ksdensity で評価する関数タイプを定義できます。または平滑化ウィンドウの帯域幅を指定できます。

[f,xi,bw] = ksdensity(___) は、カーネル平滑化ウィンドウの帯域幅 bw も返します。既定の帯域幅は、正規分布の密度に最適です。

ksdensity(___) は、カーネル平滑化関数推定をプロットします。

ksdensity(ax,___) は、gca で返される現在の Axes ではなく、ハンドル ax をもつ Axes を使用して結果をプロットします。

すべて折りたたむ

2 つの正規分布の組み合わせから標本データ セットを生成します。

rng('default')  % For reproducibility
x = [randn(30,1); 5+randn(30,1)];

推定密度をプロットします。

[f,xi] = ksdensity(x); 
figure
plot(xi,f);

Figure contains an axes object. The axes object contains an object of type line.

密度推定は、標本の二峰性を示します。

半正規分布から非負の標本データ セットを生成します。

rng('default') % For reproducibility
pd = makedist('HalfNormal','mu',0,'sigma',1);
x = random(pd,100,1);

名前と値のペアの引数 'BoundaryCorrection' を使用して、対数変換および反転という 2 つの異なる境界補正方法で pdf を推定します。

pts = linspace(0,5,1000); % points to evaluate the estimator
[f1,xi1] = ksdensity(x,pts,'Support','positive');
[f2,xi2] = ksdensity(x,pts,'Support','positive','BoundaryCorrection','reflection');

推定した 2 つの pdf をプロットします。

plot(xi1,f1,xi2,f2)
lgd = legend('log','reflection');
title(lgd, 'Boundary Correction Method')
xl = xlim;
xlim([xl(1)-0.25 xl(2)])

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent log, reflection.

ksdensity は、正または有界のいずれかのサポートが指定された場合に境界補正方法を使用します。既定の境界補正方法は対数変換です。ksdensity は、サポートを元の状態に変換するときに、カーネル密度推定量に 1/x 項を導入します。したがって、推定では x = 0 の近くにピークがあります。これに対して反転法では、望ましくないピークが境界の付近に発生することはありません。

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

load hospital

指定された一連の値で評価された推定累積分布関数を計算およびプロットします。

pts = (min(hospital.Weight):2:max(hospital.Weight));
figure()
ecdf(hospital.Weight)
hold on
[f,xi,bw] = ksdensity(hospital.Weight,pts,'Support','positive',...
	'Function','cdf');
plot(xi,f,'-g','LineWidth',2)
legend('empirical cdf','kernel-bw:default','Location','northwest')
xlabel('Patient weights')
ylabel('Estimated cdf')

Figure contains an axes object. The axes object with xlabel Patient weights, ylabel Estimated cdf contains 2 objects of type stair, line. These objects represent empirical cdf, kernel-bw:default.

ksdensity は、累積分布関数推定の平滑化が過剰であるように見えます。帯域幅を狭くして推定にすると、経験的累積分布関数に近い推定が生成される可能性があります。

平滑化ウィンドウの帯域幅を返します。

bw
bw = 0.1070

より狭い帯域幅を使用して累積分布関数推定をプロットします。

[f,xi] = ksdensity(hospital.Weight,pts,'Support','positive',...
	'Function','cdf','Bandwidth',0.05); 
plot(xi,f,'--r','LineWidth',2)
legend('empirical cdf','kernel-bw:default','kernel-bw:0.05',...
	'Location','northwest')
hold off

Figure contains an axes object. The axes object with xlabel Patient weights, ylabel Estimated cdf contains 3 objects of type stair, line. These objects represent empirical cdf, kernel-bw:default, kernel-bw:0.05.

より小さい帯域幅を使用する ksdensity 推定は、経験的累積分布関数に一致する確率が高くなります。

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

load hospital

等間隔に配置された 50 個の点で評価される推定累積分布関数をプロットします。

figure()
ksdensity(hospital.Weight,'Support','positive','Function','cdf',...
'NumPoints',50)
xlabel('Patient weights')
ylabel('Estimated cdf')

Figure contains an axes object. The axes object with xlabel Patient weights, ylabel Estimated cdf contains an object of type line.

平均値 3 をもつ指数分布から標本データを生成します。

rng('default')  % For reproducibility
x = random('exp',3,100,1);

打ち切りを示す logical ベクトルを作成します。ここで、10 よりも寿命の長い観測値は打ち切られます。

T = 10;
cens = (x>T);

推定密度関数を計算およびプロットします。

figure
ksdensity(x,'Support','positive','Censoring',cens);

Figure contains an axes object. The axes object contains an object of type line.

生存時間関数を計算およびプロットします。

figure
ksdensity(x,'Support','positive','Censoring',cens,...
'Function','survivor');

Figure contains an axes object. The axes object contains an object of type line.

累積ハザード関数を計算およびプロットします。

figure
ksdensity(x,'Support','positive','Censoring',cens,...
'Function','cumhazard');

Figure contains an axes object. The axes object contains an object of type line.

2 つの正規分布の混合を生成し、推定逆累積分布関数を指定された一連の確率値でプロットします。

rng('default')  % For reproducibility
x = [randn(30,1); 5+randn(30,1)];
pi = linspace(.01,.99,99);
figure
ksdensity(x,pi,'Function','icdf');

Figure contains an axes object. The axes object contains an object of type line.

2 つの正規分布の混合を生成します。

rng('default')  % For reproducibility
x = [randn(30,1); 5+randn(30,1)];

確率密度推定の平滑化ウィンドウの帯域幅を返します。

[f,xi,bw] = ksdensity(x); 
bw
bw = 1.5141

既定の帯域幅は、正規分布の密度に最適です。

推定密度をプロットします。

figure
plot(xi,f);
xlabel('xi')
ylabel('f')
hold on

Figure contains an axes object. The axes object with xlabel xi, ylabel f contains an object of type line.

増大した帯域幅値を使用して密度をプロットします。

[f,xi] = ksdensity(x,'Bandwidth',1.8);
plot(xi,f,'--r','LineWidth',1.5)

Figure contains an axes object. The axes object with xlabel xi, ylabel f contains 2 objects of type line.

帯域幅が広がると密度推定がさらに平滑化されるため、分布の一部の特徴が認識されない可能性があります。

次に、減少した帯域幅値を使用して密度をプロットします。

[f,xi] = ksdensity(x,'Bandwidth',0.8);
plot(xi,f,'-.k','LineWidth',1.5)
legend('bw = default','bw = 1.8','bw = 0.8')
hold off

Figure contains an axes object. The axes object with xlabel xi, ylabel f contains 3 objects of type line. These objects represent bw = default, bw = 1.8, bw = 0.8.

帯域幅が小さくなると密度推定の平滑化は少なくなり、標本のいくつかの特徴が誇張されます。

密度を評価する点が格納されている 2 列のベクトルを作成します。

gridx1 = -0.25:.05:1.25;
gridx2 = 0:.1:15;
[x1,x2] = meshgrid(gridx1, gridx2);
x1 = x1(:);
x2 = x2(:);
xi = [x1 x2];

複数の二変量正規分布の混合から取得した乱数が格納されている 30 行 2 列の行列を生成します。

rng('default')  % For reproducibility
x = [0+.5*rand(20,1) 5+2.5*rand(20,1);
            .75+.25*rand(10,1) 8.75+1.25*rand(10,1)];

標本データの推定密度をプロットします。

figure
ksdensity(x,xi);

Figure contains an axes object. The axes object contains an object of type surface.

入力引数

すべて折りたたむ

ksdensityf の値を返す対象となる標本データ。列ベクトルまたは 2 列の行列を指定します。一変量データの場合は列ベクトル、二変量データの場合は 2 列の行列を使用します。

例: [f,xi] = ksdensity(x)

データ型: single | double

f を評価する点。ベクトルまたは 2 列の行列を指定します。一変量データの場合、pts は行ベクトルまたは列ベクトルにすることができます。返される出力 f の長さは、pts 内の点の個数と同じです。

例: pts = (0:1:25); ksdensity(x,pts);

データ型: single | double

ksdensity のプロット先の Figure の Axes ハンドル。ハンドルとして指定します。

たとえば、h が Figure のハンドルである場合、ksdensity はその Figure に次のようにプロットできます。

例: ksdensity(h,x)

名前と値の引数

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

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

例: 'Censoring',cens,'Kernel','triangle','NumPoints',20,'Function','cdf' は、ksdensity が三角カーネル平滑化関数を使用し、データの範囲全体で等間隔に配置された 20 個の点を評価して累積分布関数の推定を指定し、ベクトル cens の打ち切り情報を表しています。

カーネル平滑化ウィンドウの帯域幅。x に格納されている点の数の関数です。'Bandwidth' とスカラー値から構成されるコンマ区切りのペアとして指定します。標本データが二変量の場合、Bandwidth を 2 要素のベクトルにすることもできます。既定の設定は正規分布の密度の推定に最適ですが[1]、大きい値を選択してさらに平滑化したり、小さい値を選択して平滑化を少なくすることもできます。

'BoundaryCorrection' として 'log' (既定値)、'Support' として 'positive' またはベクトル [L U] を指定した場合、ksdensity は対数変換を使用して有界のデータを非有界に変換します。'Bandwidth' の値は、変換された値のスケールに対するものになります。

例: 'Bandwidth',0.8

データ型: single | double

境界補正方法。''BoundaryCorrection'' と、'log' または 'reflection' から構成されるコンマ区切りのペアとして指定します。

説明
'log'

ksdensity は、次のいずれかの変換によって有界のデータ x を非有界に変換します。その後、密度推定の後、元の有界スケールに戻します。

  • 一変量データの場合、'Support','positive' を指定すると、ksdensitylog(x) を適用します。

  • 一変量データの場合、'Support',[L U] を指定すると (LU は数値スカラーであり L < U)、ksdensitylog((x-L)/(U–x)) を適用します。

  • 二変量データの場合、ksdensity は一変量データの場合と同じ方法で x の各列を変換します。

'Bandwidth' の値と出力 bw は、変換された値のスケールに対するものになります。

'reflection'

ksdensity は、反転データを境界の近くに追加することにより有界データを補充してから、元のサポートに対応する推定値を返します。詳細については、反転法を参照してください。

ksdensity は、'Support' として 'unbounded' 以外の値が指定された場合のみ、境界補正を適用します。

例: 'BoundaryCorrection','reflection'

打ち切るエントリを指定する logical ベクトル。'Censoring' とバイナリ値のベクトルから構成されるコンマ区切りのペアとして指定します。値が 0 の場合は打ち切りがないことを示し、1 は観測が打ち切られることを示します。既定の設定では、打ち切りはありません。この名前と値のペアの引数は、一変量データの場合のみ有効です。

例: 'Censoring',censdata

データ型: logical

推定する関数。'Function' と、以下のいずれかで構成されるコンマ区切りペアとして指定します。

説明
'pdf'確率密度関数。
'cdf'累積分布関数。
'icdf'

逆累積分布関数。ksdensityx の値について推定した逆累積分布関数を計算し、pi で指定された確率値で評価します。

この値は一変量データのみについて有効です。

'survivor'生存時間関数。
'cumhazard'

累積ハザード関数。

この値は一変量データのみについて有効です。

例: 'Function','icdf'

カーネル平滑化のタイプ。'Kernel' と、以下のいずれかで構成されるコンマ区切りペアとして指定します。

  • 'normal' (既定の設定)

  • 'box'

  • 'triangle'

  • 'epanechnikov'

  • カスタム関数または組み込み関数であるカーネル関数。関数ハンドル (@myfunction@normpdf など) を使用するか、文字ベクトルまたは string スカラー ('myfunction''normpdf' など) を使用して、関数を指定します。指定した関数は、1 つの引数を使用して呼び出されます。この引数は、密度の評価位置とデータ値の間の距離の配列です。関数は、カーネル関数の対応する値を含む、同じサイズの配列を返さなければなりません。

    'Function''pdf' である場合、カーネル関数は密度値を返します。それ以外の場合は、累積確率値を返します。

    'Function' 値が 'icdf' の場合にカスタム カーネルを指定すると、エラーが返されます。

二変量データの場合、ksdensity は同じカーネルを各次元に適用します。

例: 'Kernel','box'

xi 内の等間隔に配置された点の数。'NumPoints' とスカラー値で構成されるコンマ区切りのペアとして指定します。この名前と値のペアの引数は、一変量データの場合のみ有効です。

たとえば、標本データの範囲内で等間隔に配置された 80 個の点で、指定された関数のカーネル平滑化推定値を求めるには、次のように入力します。

例: 'NumPoints',80

データ型: single | double

密度のサポート。'support' と、以下のいずれかで構成されるコンマ区切りペアとして指定します。

説明
'unbounded'既定の設定。密度は、実数直線全体に広がります。
'positive'密度を正の値に制限します。
2 要素ベクトル [L U]密度をサポートするための下限と上限を指定します。このオプションは、一変量標本データの場合のみ有効です。
2 行 2 列の行列 [L1 L2; U1 U2]密度をサポートするための下限と上限を指定します。1 行目には下限が、2 行目には上限が格納されます。このオプションは、二変量標本データの場合のみ有効です。

二変量データの場合、'Support' には [0 -Inf; Inf Inf][0 L; Inf U] のような、正、非有界または有界の変数の組み合わせを指定できます。

例: 'Support','positive'

例: 'Support',[0 10]

データ型: single | double | char | string

カーネル密度プロットの作成に使用する関数。'PlotFcn' と次のいずれかから構成されるコンマ区切りのペアとして指定します。

説明
'surf'surf を使用して作成される 3 次元影付き表面プロット
'contour'contour を使用して作成される等高線図
'plot3'plot3 を使用して作成される 3 次元ライン プロット
'surfc'surfc を使用して作成される 3 次元影付き表面プロットの下の等高線図

この名前と値のペアの引数は、二変量標本データの場合のみ有効です。

例: 'PlotFcn','contour'

標本データの重み。'Weights' と長さが size(x,1) のベクトルから構成されるコンマ区切りのペアとして指定します。x は標本データです。

例: 'Weights',xw

データ型: single | double

出力引数

すべて折りたたむ

推定される関数値。xi または pts 内の点の個数と等しい長さをもつベクトルとして返されます。

ksdensityf を計算する評価点。ベクトルまたは 2 列の行列として返されます。一変量データの場合、既定の設定は x 内のデータ範囲全体で等間隔に配置された 100 個の点です。二変量データの場合、既定の設定は等間隔に配置された 900 個の点です。これらの点は、meshgrid を使用して各次元で等間隔に配置された 30 個の点から作成されます。

平滑化ウィンドウの帯域幅。スカラー値として返されます。

'BoundaryCorrection' として 'log' (既定値)、'Support' として 'positive' またはベクトル [L U] を指定した場合、ksdensity は対数変換を使用して有界のデータを非有界に変換します。bw の値は、変換された値のスケールに対するものになります。

詳細

すべて折りたたむ

カーネル分布

カーネル分布は、確率変数の確率密度関数 (pdf) のノンパラメトリック表現です。パラメトリック分布ではデータを適切に記述できなかったり、データの分布に関する仮定を実行しない場合、カーネル分布を使用できます。カーネル分布は、平滑化関数と帯域幅の値によって定義されます。これらは、生成される密度曲線の滑らかさを制御します。

カーネル密度推定量は、確率変数の推定 pdf です。x の任意の実数値について、カーネル密度推定量は次の式によって与えられます。

f^h(x)=1nhi=1nK(xxih),

ここで、x1、x2、...、xn は未知の分布に由来する無作為標本、n は標本サイズ、K(·) はカーネル平滑化関数、h は帯域幅です。

x の任意の実数値について、累積分布関数 (cdf) のカーネル推定量は次の式によって与えられます。

F^h(x)=xf^h(t)dt=1ni=1nG(xxih),

ここで、G(x)=xK(t)dt です。

詳細は、カーネル分布を参照してください。

反転法

反転法は、確率変数に有界のサポートがある場合にカーネル密度推定量を正確に求める境界補正方法です。'BoundaryCorrection','reflection' が指定された場合、ksdensity は反転法を使用します。この方法では、反転データを境界の近くに追加することにより有界データを補充し、pdf を推定します。そして、ksdensity は、元のサポートに対して推定 pdf を積分すると 1 になるよう適切に正規化された、元のサポートに対応する推定 pdf を返します。

さらに 'Support',[L U] が指定された場合、ksdensity はカーネル推定量を以下のように求めます。

  • 'Function''pdf' である場合、カーネル密度推定量は

    f^h(x)=1nhi=1n[K(xxih)+K(xxih)+K(xxi+h)] (L ≤ x ≤ U)。

    ここで、xi=2Lxi および xi+=2Uxi であり、xii 番目の標本データです。

  • 'Function''cdf' である場合、cdf のカーネル推定量は

    F^h(x)=1ni=1n[G(xxih)+G(xxih)+G(xxi+h)]1ni=1n[G(Lxih)+G(Lxih)+G(Lxi+h)] (L ≤ x ≤ U)。

  • ('Function''icdf''survivor'、または 'cumhazrd' である場合に) 累積分布逆関数、生存時間関数または累積ハザード関数のカーネル推定量を取得するために、ksdensityf^h(x)F^h(x) の両方を使用します。

さらに 'Support' として 'positive' または [0 inf] が指定された場合、ksdensity は上記の方程式で [L U][0 inf] に置き換えることによりカーネル推定量を求めます。

代替機能

一変量データの pdf または cdf は、MATLAB® の関数 kde を使用して推定することもできます。ksdensity と異なり、kde では境界補正方法やデータの打ち切りはサポートされません。

参照

[1] Bowman, A. W., and A. Azzalini. Applied Smoothing Techniques for Data Analysis. New York: Oxford University Press Inc., 1997.

[2] Hill, P. D. “Kernel estimation of a distribution function.” Communications in Statistics - Theory and Methods. Vol 14, Issue. 3, 1985, pp. 605-620.

[3] Jones, M. C. “Simple boundary correction for kernel density estimation.” Statistics and Computing. Vol. 3, Issue 3, 1993, pp. 135-146.

[4] Silverman, B. W. Density Estimation for Statistics and Data Analysis. Chapman & Hall/CRC, 1986.

拡張機能

バージョン履歴

R2006a より前に導入