Main Content

このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。

fitgmdist

データへの混合ガウス モデルの当てはめ

ページ内をすべて折りたたむ

構文

GMModel = fitgmdist(X,k)
GMModel = fitgmdist(X,k,Name,Value)

説明

例

GMModel = fitgmdist(X,k) は、データ (X) に近似された k 個の成分をもつ混合ガウス分布モデル (GMModel) を返します。

例

GMModel = fitgmdist(X,k,Name,Value) は、1 つ以上の Name,Value のペアの引数によって指定された追加オプションを使用して混合ガウス分布モデルを返します。

たとえば、正則化の値や共分散のタイプを指定できます。

例

すべて折りたたむ

ライブ スクリプトを開く

2 つの二変量ガウス分布の混合からデータを作成します。

mu1 = [1 2];
Sigma1 = [2 0; 0 0.5];
mu2 = [-3 -5];
Sigma2 = [1 0;0 1];
rng(1); % For reproducibility
X = [mvnrnd(mu1,Sigma1,1000); mvnrnd(mu2,Sigma2,1000)];

混合ガウス モデルを近似します。成分数を 2 に指定します。

GMModel = fitgmdist(X,2);

近似された混合ガウス モデルの等高線上にデータをプロットします。

figure
y = [zeros(1000,1);ones(1000,1)];
h = gscatter(X(:,1),X(:,2),y);
hold on
gmPDF = @(x,y) arrayfun(@(x0,y0) pdf(GMModel,[x0 y0]),x,y);
g = gca;
fcontour(gmPDF,[g.XLim g.YLim])
title('{\bf Scatter Plot and Fitted Gaussian Mixture Contours}')
legend(h,'Model 0','Model1')
hold off

Figure contains an axes object. The axes object with title blank S c a t t e r blank P l o t blank a n d blank F i t t e d blank G a u s s i a n blank M i x t u r e blank C o n t o u r s contains 3 objects of type line, functioncontour. These objects represent Model 0, Model1.

ライブ スクリプトを開く

2 つの二変量ガウス分布の混合からデータを作成します。3 番目の予測子 (最初と 2 番目の予測子の合計) を作成します。

mu1 = [1 2];
Sigma1 = [1 0; 0 1];
mu2 = [3 4];
Sigma2 = [0.5 0; 0 0.5];
rng(3); % For reproducibility
X1 = [mvnrnd(mu1,Sigma1,100);mvnrnd(mu2,Sigma2,100)];
X = [X1,X1(:,1)+X1(:,2)];

X の列は線形依存しています。このため、悪条件の共分散行列が発生することがあります。

混合ガウス モデルをデータに近似します。try / catch ステートメントを使用してエラー メッセージを管理することができます。

rng(1); % Reset seed for common start values
try
    GMModel = fitgmdist(X,2)
catch exception
    disp('There was an error fitting the Gaussian mixture model')
    error = exception.message
end
There was an error fitting the Gaussian mixture model

error = 
'Ill-conditioned covariance created at iteration 2.'

共分散推定は悪条件になっています。この結果、最適化は停止し、エラーが表示されます。

混合ガウス モデルを再度近似します。今回は正則化を使用します。

rng(3); % Reset seed for common start values
GMModel = fitgmdist(X,2,'RegularizationValue',0.1)
GMModel = 

Gaussian mixture distribution with 2 components in 3 dimensions
Component 1:
Mixing proportion: 0.536725
Mean:    2.8831    3.9506    6.8338

Component 2:
Mixing proportion: 0.463275
Mean:    0.8813    1.9758    2.8571

この場合、正則化によりアルゴリズムは解に収束します。

ライブ スクリプトを開く

混合ガウスモデルでは、データに近似する前に成分数を指定する必要があります。多くの応用で、適切な成分数を把握することが難しい可能性があります。この例では、データを調べ、主成分分析により成分数の初期推定を試行する方法を説明します。

フィッシャーのアヤメのデータセットを読み込みます。

load fisheriris
classes = unique(species)
classes = 3x1 cell
    {'setosa'    }
    {'versicolor'}
    {'virginica' }

このデータセットには、アヤメの種類のクラスが 3 つあります。分析では、この情報は不明であるとみなされます。

主成分分析を使用して、可視化するデータの次元を 2 次元に削減します。

[~,score] = pca(meas,'NumComponents',2);

1、2 または 3 個の成分を指定して、3 つの混合ガウス モデルを近似します。最適化の反復数を 1000 に増やします。ドット表記を使用して、パラメーターの最終推定を格納します。既定では、成分ごとに完全な個々の共分散を近似します。

GMModels = cell(3,1); % Preallocation
options = statset('MaxIter',1000);
rng(1); % For reproducibility

for j = 1:3
    GMModels{j} = fitgmdist(score,j,'Options',options);
    fprintf('\n GM Mean for %i Component(s)\n',j)
    Mu = GMModels{j}.mu
end
 GM Mean for 1 Component(s)

Mu = 1×2
10-14 ×

    0.1066   -0.0414


 GM Mean for 2 Component(s)

Mu = 2×2

    1.3212   -0.0954
   -2.6424    0.1909


 GM Mean for 3 Component(s)

Mu = 3×2

    0.4856   -0.1287
    1.4484   -0.0904
   -2.6424    0.1909

GMModels は 3 つの近似された gmdistribution モデルが格納されている cell 配列です。3 つの成分モデルの平均は異なります。このため、モデルはアヤメの 3 つの種類を区別していると考えられます。

近似された混合ガウス モデルの等高線上にスコアをプロットします。データセットにはラベルが含まれるため、gscatter を使用して真の数の成分を区別します。

figure
for j = 1:3
    subplot(2,2,j)
    h1 = gscatter(score(:,1),score(:,2),species);
    h = gca;
    hold on
    gmPDF = @(x,y) arrayfun(@(x0,y0) pdf(GMModels{j},[x0 y0]),x,y);
    fcontour(gmPDF,[h.XLim h.YLim],'MeshDensity',100)
    title(sprintf('GM Model - %i Component(s)',j));
    xlabel('1st principal component');
    ylabel('2nd principal component');
    if(j ~= 3)
        legend off;
    end
    hold off
end
g = legend(h1);
g.Position = [0.7 0.25 0.1 0.1];

Figure contains 3 axes objects. Axes object 1 with title GM Model - 1 Component(s) contains 4 objects of type line, functioncontour. These objects represent setosa, versicolor, virginica. Axes object 2 with title GM Model - 2 Component(s) contains 4 objects of type line, functioncontour. These objects represent setosa, versicolor, virginica. Axes object 3 with title GM Model - 3 Component(s) contains 4 objects of type line, functioncontour. These objects represent setosa, versicolor, virginica.

この 3 つの成分がある混合ガウス モデルを PCA と併用すると、アヤメの 3 つの種類が区別されるようです。

この他のオプションも、混合ガウス モデルの適切な成分数を選択するのに役立ちます。たとえば、以下のようにします。

  • 情報条件 (AIC や BIC など) を使用して、複数のモデルをさまざまな成分数で比較する。

  • evalclusters を使用してクラスター数を推定する。この方法では、Calinski-Harabasz 基準とギャップ統計などの基準がサポートされます。

ライブ スクリプトを開く

混合ガウスモデルでは、データに近似する前に成分数を指定する必要があります。多くの応用で、適切な成分数を把握することが難しい可能性があります。次の例では、AIC 統計値の近似を使用し、成分数の変化に応じて最適な混合ガウス モデルの当てはめを選択できるようにします。

2 つの二変量ガウス分布の混合からデータを作成します。

mu1 = [1 1];
Sigma1 = [0.5 0; 0 0.5];
mu2 = [2 4];
Sigma2 = [0.2 0; 0 0.2];
rng(1);
X = [mvnrnd(mu1,Sigma1,1000);mvnrnd(mu2,Sigma2,1000)];

plot(X(:,1),X(:,2),'ko')
title('Scatter Plot')
xlim([min(X(:)) max(X(:))]) % Make axes have the same scale
ylim([min(X(:)) max(X(:))])

Figure contains an axes object. The axes object with title Scatter Plot contains an object of type line.

基になるパラメーター値が不明だと仮定し、この散布図は次のことを示唆します。

  • 成分は 2 つあります。

  • クラスター間の分散は異なります。

  • クラスター内の分散は同じです。

  • クラスター内の共分散はありません。

二成分混合ガウス モデルを当てはめます。散布図の検証に基づき、共分散行列が対角となるよう指定します。statset 構造体を名前と値のペアの引数 Options の値として渡すと、最後の反復および対数尤度統計値がコマンド ウィンドウに出力されます。

options = statset('Display','final');
GMModel = fitgmdist(X,2,'CovarianceType','diagonal','Options',options);
11 iterations, log-likelihood = -4787.38

GMModel は近似された gmdistribution モデルです。

さまざまな成分数に対して、AIC を調べます。

AIC = zeros(1,4);
GMModels = cell(1,4);
options = statset('MaxIter',500);
for k = 1:4
    GMModels{k} = fitgmdist(X,k,'Options',options,'CovarianceType','diagonal');
    AIC(k)= GMModels{k}.AIC;
end

[minAIC,numComponents] = min(AIC);
numComponents
numComponents = 2

BestModel = GMModels{numComponents}
BestModel = 

Gaussian mixture distribution with 2 components in 2 dimensions
Component 1:
Mixing proportion: 0.501719
Mean:    1.9824    4.0013

Component 2:
Mixing proportion: 0.498281
Mean:    0.9880    1.0511

AIC が最も小さくなるのは、成分数が 2 の混合ガウス モデルを近似する場合です。

ライブ スクリプトを開く

混合ガウス モデルのパラメーターの推定値は、さまざまな初期値によって変化する場合があります。この例では、fitgmdist を使用して混合ガウス モデルを近似する際に初期値を制御する方法を説明します。

フィッシャーのアヤメのデータセットを読み込みます。花弁の長さと幅を予測子として使用します。

load fisheriris
X = meas(:,3:4);

既定の初期値を使用して、混合ガウス モデルをデータに近似します。アヤメの種類は 3 種類あるため、成分数として k = 3 を指定します。

rng(10); % For reproducibility
GMModel1 = fitgmdist(X,3);

既定の設定では、次の処理が行われます。

  1. 初期化のための k-means++ アルゴリズムを実装して、k = 3 個の初期クラスター中心を選択する。

  2. 初期共分散行列を対角行列として設定する (要素 (j, j) は X(:,j) の分散)。

  3. 初期混合比率を均一として扱う。

各観測値をラベルに結合し、混合ガウス モデルを近似します。

y = ones(size(X,1),1);
y(strcmp(species,'setosa')) = 2;
y(strcmp(species,'virginica')) = 3;

GMModel2 = fitgmdist(X,3,'Start',y);

初期平均、共分散行列、混合比率を明確に指定して、混合ガウス モデルを近似します。

Mu = [1 1; 2 2; 3 3];
Sigma(:,:,1) = [1 1; 1 2];
Sigma(:,:,2) = 2*[1 1; 1 2];
Sigma(:,:,3) = 3*[1 1; 1 2];
PComponents = [1/2,1/4,1/4];
S = struct('mu',Mu,'Sigma',Sigma,'ComponentProportion',PComponents);

GMModel3 = fitgmdist(X,3,'Start',S);

gscatter を使用して、アヤメの種類を区別する散布図をプロットします。モデルごとに、近似された混合ガウス モデルの等高線をプロットします。

figure
subplot(2,2,1)
h = gscatter(X(:,1),X(:,2),species,[],'o',4);
haxis = gca;
xlim = haxis.XLim;
ylim = haxis.YLim;
d = (max([xlim ylim])-min([xlim ylim]))/1000;
[X1Grid,X2Grid] = meshgrid(xlim(1):d:xlim(2),ylim(1):d:ylim(2));
hold on
contour(X1Grid,X2Grid,reshape(pdf(GMModel1,[X1Grid(:) X2Grid(:)]),...
    size(X1Grid,1),size(X1Grid,2)),20)
uistack(h,'top')
title('{\bf Random Initial Values}');
xlabel('Sepal length');
ylabel('Sepal width');
legend off;
hold off
subplot(2,2,2)
h = gscatter(X(:,1),X(:,2),species,[],'o',4);
hold on
contour(X1Grid,X2Grid,reshape(pdf(GMModel2,[X1Grid(:) X2Grid(:)]),...
    size(X1Grid,1),size(X1Grid,2)),20)
uistack(h,'top')
title('{\bf Initial Values from Labels}');
xlabel('Sepal length');
ylabel('Sepal width');
legend off
hold off
subplot(2,2,3)
h = gscatter(X(:,1),X(:,2),species,[],'o',4);
hold on
contour(X1Grid,X2Grid,reshape(pdf(GMModel3,[X1Grid(:) X2Grid(:)]),...
    size(X1Grid,1),size(X1Grid,2)),20)
uistack(h,'top')
title('{\bf Initial Values from the Structure}');
xlabel('Sepal length');
ylabel('Sepal width');
legend('Location',[0.7,0.25,0.1,0.1]);
hold off

Figure contains 3 axes objects. Axes object 1 with title blank R a n d o m blank I n i t i a l blank V a l u e s contains 4 objects of type contour, line. These objects represent virginica, versicolor, setosa. Axes object 2 with title blank I n i t i a l blank V a l u e s blank f r o m blank L a b e l s contains 4 objects of type contour, line. These objects represent virginica, versicolor, setosa. Axes object 3 with title blank I n i t i a l blank V a l u e s blank f r o m blank t h e blank S t r u c t u r e contains 4 objects of type contour, line. These objects represent virginica, versicolor, setosa.

等高線から、GMModel2 が軽い三峰性を、他が二峰性分布を示していることがわかります。

推定成分平均値を表示します。

table(GMModel1.mu,GMModel2.mu,GMModel3.mu,'VariableNames',...
    {'Model1','Model2','Model3'})
ans=3×3 table
         Model1               Model2              Model3     
    _________________    ________________    ________________

    5.2115     2.0119    4.2857    1.3339    1.4604    0.2429
     1.461    0.24423     1.462     0.246    4.7509    1.4629
    4.6829     1.4429    5.5507    2.0316    5.0158    1.8592

アヤメの種類を最もよく区別しているのは GMModel2 であると判断できます。

入力引数

すべて折りたたむ

混合ガウス モデルを近似させるデータ。数値行列として指定します。

X の行は観測値に、X の列は変数に対応します。観測値の個数は、変数の個数および成分の個数より多くなければなりません。

NaN は欠損値を表します。近似の前に、X の行から少なくとも 1 つの NaN をもつ行が削除されます。これにより、有効な標本サイズが小さくなります。

データ型: single | double

混合ガウス モデルを近似させる際に使用する成分の数。正の整数として指定します。たとえば、k = 3 を指定すると、3 つの異なる平均値、共分散行列、成分比率をもつ混合ガウス モデルでデータ (X) が近似されます。

データ型: single | double

名前と値の引数

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

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

例: 'RegularizationValue',0.1,'CovarianceType','diagonal' は、正則化パラメーターの値を 0.1 に設定し、対角共分散行列を近似するよう指定します。

データに近似させる分散共分散行列のタイプ。'CovarianceType' と、'diagonal' または 'full' のいずれかで構成される、コンマ区切りのペアとして指定します。

'diagonal' を設定すると、対角共分散行列が近似されます。この場合、k*d 個の共分散パラメーターが推定されます (d は X の列数、すなわち d = size(X,2))。

それ以外の場合、完全な共分散行列が近似されます。この場合、k*d*(d+1)/2 個の共分散パラメーターが推定されます。

例: 'CovarianceType','diagonal'

反復 EM アルゴリズム最適化オプション。'Options' と statsetオプション構造体で構成されるコンマ区切りのペアとして指定します。

使用可能な名前と値のペアの引数を以下の表に示します。

名前値
'Display'

'final':最終出力を表示します。

'iter':一部の関数ではコマンド ウィンドウに反復出力が表示されます。それ以外では最終出力が表示されます。

'off':最適化情報は表示されません。

'MaxIter'許可される反復の最大数を示す正の整数。既定値は 100 です。
'TolFun'対数尤度関数値の終了許容誤差を示す正の整数。既定の設定は 1e-6 です。

例: 'Options',statset('Display','final','MaxIter',1500,'TolFun',1e-5)

事後確率の許容誤差。ProbabilityTolerance と範囲 [0,1e-6] の非負のスカラー値から構成されるコンマ区切りのペアとして指定します。

各反復で事後確率を推定した後で、fitgmdist は許容誤差値以下の事後確率をゼロに設定します。非ゼロの許容誤差を使用すると、fitgmdist が高速になる可能性があります。

例: 'ProbabilityTolerance',0.0000025

データ型: single | double

正則化パラメーター値。'RegularizationValue' と非負のスカラーで構成されるコンマ区切りのペアとして指定します。

推定される共分散行列を正定値行列にするには、RegularizationValue を小さい正のスカラーに設定します。

例: 'RegularizationValue',0.01

データ型: single | double

初期値の新しいセットを使用して EM アルゴリズムを繰り返す回数。'Replicates' と正の整数で構成されるコンマ区切りのペアとして指定します。

Replicates が 1 より大きい場合は、次のようになります。

  • 名前と値のペアの引数 Start は plus (既定値) または randSample でなければなりません。

  • GMModel は最も大きい対数尤度で近似されます。

例: 'Replicates',10

データ型: single | double

すべての共分散行列が同一である (プールされた推定を近似する) かどうかを示すフラグ。'SharedCovariance' と、論理値 false または true で構成される、コンマ区切りのペアとして指定します。

SharedCovariance が true の場合、k 個の共分散行列はすべて等しく、共分散パラメーターの数は k 分の 1 に減ります。

初期値の設定方法。'Start' と 'randSample'、'plus'、整数のベクトルまたは構造体配列で構成されるコンマ区切りのペアとして指定します。

Start の値により、各ガウス成分パラメーター (平均、共分散、混合比率) の最適化ルーチンで必要となる初期値が決定されます。次の表は、使用できるオプションの一覧です。

値説明
'randSample'X から k 個の観測値が、初期成分の平均として無作為に選択されます。混在の比率は統一されています。すべてのコンポーネントの初期共分散行列は対角行列です。対角行列の要素 j は、X(:,j) の分散です。
'plus'k-means++ アルゴリズムを使用して X から k の観測値が選択されます。初期値の混在の比率は統一されています。すべてのコンポーネントの初期共分散行列は対角行列です。対角行列の要素 j は、X(:,j) の分散です。
整数のベクトル 各点の成分インデックスの初期推定が格納された、長さが n (観測値の数) のベクトル。つまり、各要素は 1 ～ k の範囲の整数であり、成分と一致しています。同じ成分に対応するすべての観測値が収集され、それぞれの平均、共分散、混合比率が計算され、初期値がこれらの統計値に設定されます。
構造体配列

d 個の変数があるとします (d = size(X,2))。構造体配列 (S など) には次の 3 つのフィールドがなければなりません。

  • S.mu:各成分の初期平均を指定する k 行 d 列の行列

  • S.Sigma:各成分の共分散行列を指定する数値配列。Sigma は次のいずれかです。

    • d x d x k の配列。Sigma(:,:,j) は成分 j の初期共分散行列です。

    • 1 x d x k の配列。diag(Sigma(:,:,j)) は成分 j の初期共分散行列です。

    • d 行 d 列の行列。Sigma はすべての成分の初期共分散行列です。

    • 1 行 d 列のベクトル。diag(Sigma) はすべての成分の初期共分散行列です。

  • S.ComponentProportion:各成分の初期混合比率を指定する 1 行 k 列のスカラーのベクトル。既定値は uniform です。

例: 'Start',ones(n,1)

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

出力引数

すべて折りたたむ

近似された混合ガウス モデル。gmdistribution モデルとして返します。

GMModel のプロパティにはドット表記でアクセスします。たとえば、GMModel.AIC と入力すると AIC が表示されます。

ヒント

fitgmdist では次の処理が行われる場合があります。

  • 1 つ以上の成分が、悪条件または特異共分散行列をもつような解に収束することがあります。

    次の問題があると、悪条件の共分散行列になる可能性があります。

    • データの次元数が比較的多く、十分な観測値がない場合。

    • データの一部の予測子 (変数) に高い相関関係がある場合。

    • 特徴量の一部またはすべてが離散である場合。

    • データの近似対象となる成分が多すぎる場合。

    通常、次のいずれかの対処法により、悪条件の共分散行列の発生を回避できます。

    • データを事前処理し、相関関係のある特徴量を排除します。

    • 'SharedCovariance' を true に設定し、すべての成分に同じ共分散行列を使用します。

    • 'CovarianceType' を 'diagonal' に設定します。

    • 'RegularizationValue' を使用し、共分散行列ごとの対角行列に対して非常に小さい正の数値を追加します。

    • 別のセットの初期値を試します。

  • 1 つ以上の成分が、悪条件の共分散行列をもつ中間ステップを経由する。別の初期値のセットを試し、データまたはモデルを変更せずにこの問題を回避してください。

アルゴリズム

すべて折りたたむ

混合ガウス モデルの尤度の最適化

混合ガウス モデルの尤度の最適化には、反復期待値最大化 (EM) アルゴリズムを使用します。

fitgmdist では、反復的な "期待値最大化" (EM) アルゴリズムを使用して GMM をデータに当てはめます。EM アルゴリズムでは、成分の平均、共分散行列および混合比率の初期値を使用し、次のステップに従って処理を進めます。

  1. 各観測値について、成分メンバーシップの事後確率を計算します。この結果は、観測値 i が成分 j から派生する事後確率が要素 (i,j) に格納されている n 行 k 列の行列と考えることができます。これは、EM アルゴリズムの "E" のステップです。

  2. 成分メンバーシップの事後確率を重みとして使用し、最大尤度を適用することにより成分の平均、共分散行列および混合比率を推定します。これは、EM アルゴリズムの "M" のステップです。

これらのステップは、収束するまで繰り返されます。尤度面は複雑なので、局所的な最適解に収束する場合があります。また、結果の局所的な最適解は初期条件によって異なる場合があります。fitgmdist には、k-means++ アルゴリズムや観測値についての無作為な成分割り当てなど、初期条件を選択するためのオプションがいくつかあります。

初期化のための k-means++ アルゴリズム

k-means++ アルゴリズムはヒューリスティックを使用して k-means クラスタリングの重心シードを検出します。fitgmdist でも同じ原則を適用し、k-means++ アルゴリズムを使用して EM アルゴリズムを初期化し、近似混合ガウス モデルの初期パラメーター値を選択できます。

k-means++ アルゴリズムはクラスター数を k であると仮定し、初期パラメーター値を以下のように選択します。

  1. 成分混合確率として、次の一様確率を選択します。pi=1kここで、i = 1、...、k です。

  2. 共分散行列として、同一の対角行列を選択します。ここで、σi=diag(a1,a2,…,ak) および aj=var(Xj) です。

  3. 最初の初期成分の中心 μ1 を X のすべてのデータ点から等間隔で選択します。

  4. 中心 j を選択するには、以下を行います。

    1. 各観測から各重心までのマハラノビス距離を計算します。各観測を最も近い重心に割り当てます。

    2. m = 1,...,n および p = 1,...,j - 1 について、次の確率で X から重心 j を無作為に選択します。

      d2(xm,μp)∑h;xh∈Μpd2(xh,μp)

      ここで、d(xm,μp) は観測値 m と μp の間の距離、Mp は重心 μp に最も近いすべての観測値の集合です。xm は Mp に属します。

      つまり、それ自体から、既に選択されたなかで最も近い中心までの距離に比例する確率で、次の中心を選択します。

  5. k 個の重心が選択されるまでステップ 4 を繰り返します。

参照

[1] McLachlan, G., and D. Peel. Finite Mixture Models. Hoboken, NJ: John Wiley & Sons, Inc., 2000.

バージョン履歴

R2014a で導入

参考

|

トピック