fitgmdist
データへの混合ガウス モデルの当てはめ
説明
例
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
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 = 3×1 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-15 ×
0.9607 -0.4371
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];
この 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(:))])
基になるパラメーター値が不明だと仮定し、この散布図は次のことを示唆します。
成分は 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);既定の設定では、次の処理が行われます。
初期化のための k-means++ アルゴリズムを実装して、k = 3 個の初期クラスター中心を選択する。
初期共分散行列を対角行列として設定する (要素 (
j,j) はX(:,j)の分散)。初期混合比率を均一として扱う。
各観測値をラベルに結合し、混合ガウス モデルを当てはめます。
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
等高線から、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 であると判断できます。
入力引数
名前と値の引数
出力引数
ヒント
fitgmdist では次の処理が行われる場合があります。
1 つ以上の成分が、悪条件または特異共分散行列をもつような解に収束することがあります。
次の問題があると、悪条件の共分散行列になる可能性があります。
データの次元数が比較的多く、十分な観測値がない場合。
データの一部の予測子 (変数) に高い相関関係がある場合。
特徴量の一部またはすべてが離散である場合。
データの近似対象となる成分が多すぎる場合。
通常、次のいずれかの対処法により、悪条件の共分散行列の発生を回避できます。
データを事前処理し、相関関係のある特徴量を排除します。
'SharedCovariance'をtrueに設定し、すべての成分に同じ共分散行列を使用します。'CovarianceType'を'diagonal'に設定します。'RegularizationValue'を使用し、共分散行列ごとの対角行列に対して非常に小さい正の数値を追加します。別のセットの初期値を試します。
1 つ以上の成分が、悪条件の共分散行列をもつ中間ステップを経由する。別の初期値のセットを試し、データまたはモデルを変更せずにこの問題を回避してください。
アルゴリズム
参照
[1] McLachlan, G., and D. Peel. Finite Mixture Models. Hoboken, NJ: John Wiley & Sons, Inc., 2000.
バージョン履歴
R2014a で導入
