ksdensity

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

構文

[f,xi] = ksdensity(x)

[f,xi] = ksdensity(x,pts)

[f,xi] = ksdensity(___,Name,Value)

[f,xi,bw] = ksdensity(___)

ksdensity(___)

ksdensity(ax,___)

説明

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

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

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

[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);

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

境界補正を使用した密度の推定

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

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)])

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')

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

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

指定された数の点の推定累積分布関数をプロットする

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

load hospital

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

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

打ち切られた故障データの生存時間および累積ハザードを推定する

平均値 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
ksdensity(x,'Support','positive','Censoring',cens,...
'Function','survivor');

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

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

指定された確率値で逆累積分布関数を推定する

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

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

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

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

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

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

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

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

[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

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

二変量データのカーネル密度推定値のプロット

密度を評価する点が格納されている 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);

入力引数

x — 標本データ

列ベクトル | 2 列の行列

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

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

データ型: single | double

pts — f を評価する点

ベクトル | 2 列の行列

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

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

データ型: single | double

ax — Axes ハンドル

ハンドル

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 の打ち切り情報を表しています。

Bandwidth — カーネル平滑化ウィンドウの帯域幅

"normal-approx" (既定値) | "plug-in" | スカラー値 | 2 要素ベクトル

カーネル平滑化ウィンドウの帯域幅。x に格納されている点の数の関数です。"normal-approx" または "plug-in" として指定します。標本データが二変量の場合、Bandwidth を 2 要素のベクトル [L,U] にすることもできます。

• Bandwidth が "normal-approx" の場合、ksdensity は正規近似法または "Silverman の経験則" を使用して帯域幅を計算します。計算された値は正規分布の密度の推定に最適ですが[2]、大きい値を指定して平滑化を増やしたり、小さい値を指定して減らしたりすることもできます。

• Bandwidth が "plug-in" の場合、ksdensity は[1]で説明されている改善されたプラグイン法を使用して帯域幅を計算します。このプラグイン法は "Sheather-Jones" 法と呼ばれることもあります。

• Bandwidth が正のスカラーの場合、その値で確率関数の推定の滑らかさが制御されます。値が大きくなると確率関数の推定が滑らかになります。

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

例: Bandwidth=0.8

データ型: single | double

BoundaryCorrection — 境界補正方法

'log' (既定値) | 'reflection'

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

値	説明
'log'	ksdensity は、次のいずれかの変換によって有界のデータ x を非有界に変換します。その後、密度推定の後、元の有界スケールに戻します。 一変量データの場合、'Support','positive' を指定すると、ksdensity は log(x) を適用します。 一変量データの場合、'Support',[L U] を指定すると (L と U は数値スカラーであり L < U)、ksdensity は log((x-L)/(U–x)) を適用します。 二変量データの場合、ksdensity は一変量データの場合と同じ方法で x の各列を変換します。 'Bandwidth' の値と出力 bw は、変換された値のスケールに対するものになります。
'reflection'	ksdensity は、反転データを境界の近くに追加することにより有界データを補充してから、元のサポートに対応する推定値を返します。詳細については、反転法を参照してください。

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

例: 'BoundaryCorrection','reflection'

Censoring — logical ベクトル

0 のベクトル (既定値) | 0 と 1 のベクトル

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

例: 'Censoring',censdata

データ型: logical

Function — 推定する関数

'pdf' (既定値) | 'cdf' | 'icdf' | 'survivor' | 'cumhazard'

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

値	説明
'pdf'	確率密度関数。
'cdf'	累積分布関数。
'icdf'	逆累積分布関数。ksdensity は x の値について推定した逆累積分布関数を計算し、pi で指定された確率値で評価します。 この値は一変量データのみについて有効です。
'survivor'	生存時間関数。
'cumhazard'	累積ハザード関数。 この値は一変量データのみについて有効です。

例: 'Function','icdf'

Kernel — カーネル平滑化のタイプ

'normal' (既定値) | 'box' | 'triangle' | 'epanechnikov' | 関数ハンドル | 文字ベクトル | string スカラー

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

• 'normal' (既定の設定)

• 'box'

• 'triangle'

• 'epanechnikov'

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

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

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

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

例: 'Kernel','box'

NumPoints — 等間隔に配置された点の数

100 (既定値) | スカラー値

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

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

例: 'NumPoints',80

データ型: single | double

Support — 密度のサポート

"unbounded" (既定値) | "positive" | "nonnegative" | "negative" | 2 要素ベクトル | 2 行 2 列の行列

密度のサポート。次の値のいずれかとして指定します。

値	説明
"unbounded"	密度が実数直線全体に広がることを許可します (既定の設定)。
"positive"	密度を正の値に制限します。
"nonnegative"	密度を正の値と 0 に制限します。
"negative"	密度を負の値に制限します。
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' (既定値) | 'contour' | 'plot3' | 'surfc'

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

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

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

例: 'PlotFcn','contour'

Weights — 標本データの重み

ベクトル

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

例: 'Weights',xw

データ型: single | double

出力引数

f — 推定される関数値

ベクトル

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

xi — 評価点

等間隔に配置された 100 個の点 | 等間隔に配置された 900 個の点 | ベクトル | 2 列の行列

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

bw — 平滑化ウィンドウの帯域幅

スカラー値

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

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

詳細

カーネル分布

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

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

$${\widehat{f}}_h\left(x\right)=\frac{1}{nh}{\displaystyle \sum _{i=1}^{n}K\left(\frac{x-x_i}{h}\right)},$$

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

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

$${\widehat{F}}_h\left(x\right)={\displaystyle {\int }_{-\infty }^x{\widehat{f}}_h(t)dt}=\frac{1}{n}{\displaystyle \sum _{i=1}^{n}G\left(\frac{x-x_i}{h}\right)},$$

ここで、$$G(x)={\displaystyle {\int }_{-\infty }^{x}K(t)dt}$$ です。

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

反転法

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

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

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

  $${\widehat{f}}_h(x)=\frac{1}{nh}{\displaystyle \sum _{i=1}^n\left[K\left(\frac{x-x_i^{-}}{h}\right)+K\left(\frac{x-x_i}{h}\right)+K\left(\frac{x-x_i^{+}}{h}\right)\right]}$$ (L ≤ x ≤ U)。

  ここで、$$x_i^{-}=2L-x_i$$ および $$x_i^{+}=2U-x_i$$ であり、x_i は i 番目の標本データです。

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

  $${\widehat{F}}_h(x)=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left[G\left(\frac{x-x_i^{-}}{h}\right)+G\left(\frac{x-x_i}{h}\right)+G\left(\frac{x-x_i^{+}}{h}\right)\right]}-\frac{1}{n}{\displaystyle \sum _{i=1}^n\left[G\left(\frac{L-x_i^{-}}{h}\right)+G\left(\frac{L-x_i}{h}\right)+G\left(\frac{L-x_i^{+}}{h}\right)\right]}$$ (L ≤ x ≤ U)。

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

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

代替機能

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

参照

[1] Botev, Z. I., J. F. Grotowski, and D. P. Kroese. "Kernel Density Estimation via Diffusion." The Annals of Statistics, vol. 38, no. 5 (October 1, 2010). https://projecteuclid.org/journals/annals-of-statistics/volume-38/issue-5/Kernel-density-estimation-via-diffusion/10.1214/10-AOS799.full

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

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

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

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

拡張機能

tall 配列

メモリの許容量を超えるような多数の行を含む配列を計算します。

この関数は、メモリ超過のデータに対し tall 配列をサポートしますが、いくつかの制限があります。

• 入力データの受け渡しまたは並べ替えがさらに必要となる一部のオプションはサポートされていません。

  • 'BoundaryCorrection'

  • 'Censoring'

  • 'Support' (サポートは常に非有界)。

• (中央絶対偏差ではなく) 標準偏差を使用して帯域幅を計算します。

詳細は、メモリに収まらないデータの tall 配列を参照してください。

C/C++ コード生成

MATLAB® Coder™ を使用して C および C++ コードを生成します。

使用上の注意事項および制限事項:

• プロットはサポートされません。

• 名前と値のペアの引数に含まれる名前はコンパイル時の定数でなければなりません。

• 名前と値のペアの引数 'BoundaryCorrection'、'Function' および 'Kernel' の値も、コンパイル時の定数でなければなりません。たとえば、生成されたコードで名前と値のペアの引数 'Function','cdf' を使用するには、{coder.Constant('Function'),coder.Constant('cdf')} を codegen の -args の値に含めます。

• 名前と値のペアの引数 'Kernel' の値をカスタム関数ハンドルにすることはできません。カスタム カーネル関数を指定するには、文字ベクトルまたは string スカラーを使用します。

• 名前と値のペアの引数 'Support' の値について、コンパイル時のデータ型と実行時のデータ型が一致しなければなりません。

コード生成の詳細については、コード生成の紹介および一般的なコード生成のワークフローを参照してください。

GPU 配列

Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

この関数は、GPU 配列を完全にサポートします。詳細は、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2006a より前に導入

参考

histogram | mvksdensity | kde

トピック

• ksdensity を使用したカーネル分布の当てはめ
• ksdensity を使用したグループ化データへの分布の当てはめ
• 確率分布の操作
• ノンパラメトリックな経験的確率分布
• サポートされている分布