mvksdensity

多変量データのカーネル平滑化関数推定値

構文

f = mvksdensity(x,pts,'Bandwidth',bw)

f = mvksdensity(x,pts,'Bandwidth',bw,Name,Value)

説明

f = mvksdensity(x,pts,'Bandwidth',bw) は、帯域幅の値に必要な名前と値のペアの引数 bw を使用し、n 行 d 列の行列 x に格納されている標本データの確率密度推定値を pts 内の点で評価して計算します。推定は積ガウス カーネル関数に基づきます。

一変量または二変量データの場合は、代わりに ksdensity を使用します。

f = mvksdensity(x,pts,'Bandwidth',bw,Name,Value) は、1 つ以上の Name,Value ペア引数で指定された追加オプションを使用して、前の出力引数のいずれかを返します。たとえば、確率密度、累積確率、生存時間関数など、mvksdensity が評価する関数のタイプを定義できます。入力値に重みを割り当てることもできます。

例

多変量カーネル密度の推定

Hald セメント データを読み込みます。

load hald

これは、13 種類のセメント組成に対する硬化熱を測定したデータです。予測子行列 ingredients には、4 種類のセメント原料のそれぞれについての組成率が格納されています。応答行列 heat には、180 日後の硬化熱 (cal/g) が格納されています。

ingredients の最初の 3 つの観測値についてカーネル密度を推定します。

xi = ingredients(1:3,:);
f = mvksdensity(ingredients,xi,'Bandwidth',0.8);

グリッドの使用による多変量カーネル密度の推定

Hald セメント データを読み込みます。

load hald

これは、13 種類のセメント組成に対する硬化熱を測定したデータです。予測子行列 ingredients には、4 種類のセメント原料のそれぞれについての組成率が格納されています。応答行列 heat には、180 日後の硬化熱 (cal/g) が格納されています。

密度を推定する点の配列を作成します。まず、各次元で同様の数の点を使用して、各変数の範囲と間隔を定義します。

gridx1 = 0:2:22;
gridx2 = 20:5:80;
gridx3 = 0:2:24;
gridx4 = 5:5:65;

次に、ndgrid を使用して、定義した範囲および間隔による点のフル グリッドを生成します。

[x1,x2,x3,x4] = ndgrid(gridx1,gridx2,gridx3,gridx4);

最後に、変換と連結を行って、密度を推定する点が格納されている配列を作成します。この配列には、各変数について 1 つずつの列があります。

x1 = x1(:,:)';
x2 = x2(:,:)';
x3 = x3(:,:)';
x4 = x4(:,:)';
xi = [x1(:) x2(:) x3(:) x4(:)];

密度を推定します。

f = mvksdensity(ingredients,xi,...
	'Bandwidth',[4.0579 10.7345 4.4185 11.5466],...
	'Kernel','normpdf');

xi と f のサイズを表示して、mvksdensity が xi 内の各点で密度を計算したことを確認します。

size_xi = size(xi)

size_xi = 1×2

       26364           4

size_f = size(f)

size_f = 1×2

       26364           1

入力引数

x — 標本データ

数値行列

mvksdensity が確率密度の推定値を返す対象となる標本データ。n 行 d 列の数値行列を指定します。n は x 内のデータ点 (行) の数、d は次元 (列) の数です。

データ型: single | double

pts — f を評価する点

行列

確率密度推定値 f を評価する点。x と同じ列数の行列を指定します。返される推定値 f および pts は、同じ行数になります。

データ型: single | double

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

スカラー値 | d 要素のベクトル

カーネル平滑化ウィンドウの帯域幅の値。スカラー値または d 要素のベクトルを指定します。d は標本データ x 内の次元 (列) の数です。bw がスカラー値の場合、すべての次元に適用されます。

'BoundaryCorrection' として 'log' (既定値)、'Support' として 'positive' または 2 行の行列を指定した場合、mvksdensity は対数変換を使用して有界のデータを非有界に変換します。bw の値は、変換された値のスケールに対するものになります。

帯域幅に関する Silverman の経験則は次のとおりです。

$$b_i={\sigma }_i{\left\{\frac{4}{\left(d+2\right)n}\right\}}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\left(d+4\right)$}\right.},i=1,2,\mathrm{...},d,$$

ここで、d は次元数、n は観測値数、$${\sigma }_i$$ は i 番目の変量の標準偏差です [4]。

例: 'Bandwidth',0.8

データ型: single | double

名前と値の引数

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

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

例: 'Kernel','triangle','Function,'cdf' は、mvksdensity が三角カーネル関数を使用して標本データの cdf を推定することを指定します。

BoundaryCorrection — 境界補正方法

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

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

値	説明
'log'	関数 mvksdensity は、次のいずれかの変換を使用して有界のデータを非有界に変換します。その後、密度推定の後、元の有界スケールに戻します。 'Support','positive' を指定した場合、mvksdensity は各次元について log(xj) を適用します。xj は入力引数 x の j 番目の列です。 各次元の下限と上限から成る 2 行の行列を 'Support' に指定した場合、mvksdensity は各次元について log((xj-Lj)/(Uj-xj)) を適用します。Lj と Uj はそれぞれ j 番目の次元の下限と上限です。 bw の値は、変換された値のスケールに対するものになります。
'reflection'	mvksdensity は、反転データを境界の近くに追加することにより有界データを補充してから、元のサポートに対応する推定値を返します。詳細については、反転法を参照してください。

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

例: 'BoundaryCorrection','reflection'

Function — 推定する関数

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

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

値	説明
'pdf'	確率密度関数
'cdf'	累積分布関数
'survivor'	生存時間関数

例: 'Function','cdf'

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

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

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

値	説明
'normal'	正規 (ガウス) カーネル
'box'	ボックス カーネル
'triangle'	三角カーネル
'epanechnikov'	Epanechnikov カーネル

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

mvksdensity は各次元に同じカーネルを適用します。

例: 'Kernel','box'

Support — 密度のサポート

'unbounded' (既定値) | 'positive' | 2 行 d 列の行列

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

値	説明
'unbounded'	密度が数直線全体に広がることを許可します
'positive'	密度を正の値に制限します
2 行 d 列の行列	密度のサポートについて有限の下限と上限を指定します。1 行目には下限が、2 行目には上限が格納されます。各列には、x の 1 つの次元に対する範囲が格納されます。

'Support' は、[0 -Inf L; Inf Inf U] のように指定される正、有界または非有界の値の組み合わせにすることもできます。

例: 'Support','positive'

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

Weights — 標本データの重み

ベクトル

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

例: 'Weights',xw

データ型: single | double

出力引数

f — 推定される関数値

ベクトル

推定される関数値。ベクトルとして返されます。f と pts は同じ行数になります。

詳細

多変量カーネル分布

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

多変量カーネル密度推定量は、乱数ベクトルの推定 pdf です。x = (x₁, x₂, …, x_d)' は密度関数が f である d 次元の乱数ベクトル、y_i = (y_i1, y_i2, …, y_id)' は i = 1, 2, …, n について f から抽出した無作為標本であるとします。ここで、n は無作為標本の個数です。x の任意の実数ベクトルについて、多変量カーネル密度推定量は次によって与えられます。

$${\widehat{f}}_H\left(x\right)=\frac{1}{n}{\displaystyle \sum _{i=1}^n{K}_H\left(x-y_i\right)},$$

ここで、$$K_H\left(x\right)={\left|H\right|}^{-1/2}K(H^{-1/2}x)$$, $$K(·)$$ はカーネル平滑化関数、H は d 行 d 列の帯域幅の行列です。

mvksdensity は、帯域幅の対角行列と積カーネルを使用します。つまり、H^1/2 はベクトル (h₁, h₂, …, h_d) の要素が主対角にある正方対角行列です。K(x) は K(x) = k(x₁)k(x₂) ⋯k(x_d) という積の形式になります。$$k(·)$$ は 1 次元のカーネル平滑化関数です。すると、多変量カーネル密度推定量は次のようになります。

$${\widehat{f}}_H\left(x\right)=\frac{1}{n}{\displaystyle \sum _{i=1}^n{K}_H\left(x-y_i\right)}=\frac{1}{n{h}_1{h}_2\cdots h_d}{\displaystyle \sum _{i=1}^{n}K\left(\frac{x_1-y_{i1}}{h_1},\frac{x_2-y_{i2}}{h_2},\cdots ,\frac{x_d-y_{id}}{h_d}\right)}=\frac{1}{n{h}_1{h}_2\cdots h_d}{\displaystyle \sum _{i=1}^n{\displaystyle \prod _{j=1}^{d}k\left(\frac{x_j-y_{ij}}{h_j}\right)}}.$$

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

$${\widehat{F}}_H\left(x\right)={\displaystyle {\int }_{-\infty }^{x_1}{\displaystyle {\int }_{-\infty }^{x_2}\cdots {\displaystyle {\int }_{-\infty }^{x_d}{\widehat{f}}_H(t)d{t}_d\cdots }d{t}_2}d{t}_1}=\frac{1}{n}{\displaystyle \sum _{i=1}^n{\displaystyle \prod _{j=1}^{d}G\left(\frac{x_j-y_{ij}}{h_j}\right)}},$$

ここで、$$G(x_j)={\displaystyle {\int }_{-\infty }^{x_j}k(t_j)d{t}_j}$$ です。

反転法

反転法は、確率変数に有界のサポートがある場合にカーネル密度推定量を正確に求める境界補正方法です。'BoundaryCorrection','reflection' が指定された場合、mvksdensity は反転法を使用します。

さらに 'Support' に、各次元の下限と上限から成る 2 行の行列を指定した場合、mvksdensity はカーネル推定量を次のように求めます。

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

  $${\widehat{f}}_H\left(x\right)=\frac{1}{n{h}_1{h}_2\cdots h_d}{\displaystyle \sum _{i=1}^n{\displaystyle \prod _{j=1}^d\left[k\left(\frac{x_j-y_{ij}^{-}}{h_j}\right)+k\left(\frac{x_j-y_{ij}}{h_j}\right)+k\left(\frac{x_j-y_{ij}^{+}}{h_j}\right)\right]}}$$ (L_j ≤ x_j ≤ U_j)。

  ここで、$$y_{ij}^{-}=2L_j-y_{ij}$$ および $$y_{ij}^{+}=2U_j-y_{ij}$$ であり、y_ij は入力引数 x の x(i,j) に対応する i 番目の標本データの j 番目の要素です。L_j と U_j はそれぞれ、j 番目の次元の下限と上限です。

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

  $${\widehat{F}}_H(x)=\frac{1}{n}{\displaystyle \sum _{i=1}^n{\displaystyle \prod _{j=1}^d\left[G\left(\frac{x_j-y_{ij}^{-}}{h_j}\right)+G\left(\frac{x_j-y_{ij}}{h_j}\right)+G\left(\frac{x_j-y_{ij}^{+}}{h_j}\right)-G\left(\frac{L_j-y_{ij}^{-}}{h_j}\right)-G\left(\frac{L_j-y_{ij}}{h_j}\right)-G\left(\frac{L_j-y_{ij}^{+}}{h_j}\right)\right]}}$$ (L_j ≤ x_j ≤ U_j)。

• 生存時間関数のカーネル推定量を取得するため ('Function' が 'survivor' である場合)、mvksdensity は $${\widehat{f}}_H(x)$$ と $${\widehat{F}}_H(x)$$ の両方を使用します。

さらに 'Support' に 'positive'、または [0 inf] を含む行列を指定した場合、mvksdensity は上記の方程式で [Lj Uj] を [0 inf] に置き換えることによりカーネル密度推定量を求めます。

参照

[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.

[5] Scott, D. W. Multivariate Density Estimation: Theory, Practice, and Visualization. John Wiley & Sons, 2015.

拡張機能

C/C++ コード生成

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

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

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

• 名前と値のペアの引数 '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)を参照してください。

バージョン履歴

R2016a で導入

参考

ksdensity

トピック

• 確率分布の操作
• ノンパラメトリックな経験的確率分布
• サポートされている分布