Создание данных из смеси двух двумерных гауссовых распределений.

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

Поместите модель гауссовой смеси. Укажите, что существует два компонента.

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

Создание данных из смеси двух двумерных гауссовых распределений. Создайте третий предиктор, который является суммой первого и второго предикторов.

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 3.'

Оценки ковариации плохо обусловлены. Следовательно, оптимизация прекращается, и появляется ошибка.

Снова установите модель гауссовой смеси, но используйте регуляризацию.

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

Набор данных содержит три класса видов радужки. Анализ продолжается так, как будто это неизвестно.

Используйте анализ основных компонентов, чтобы уменьшить размер данных до двух размеров для визуализации.

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

Подгоните три модели смеси Гаусса к данным, указав компоненты 1, 2 и 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.1044   -0.0429

 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 - массив ячеек, содержащий три gmdistribution модели. Средства в трех компонентных моделях различны, что позволяет предположить, что модель различает три вида радужки.

Постройте графики по границам подогнанной модели гауссовой смеси. Поскольку набор данных включает метки, используйте 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];

Трехкомпонентная модель гауссовой смеси в сочетании с РСА выглядит так, будто различает три вида радужки.

Есть и другие опции, которые можно использовать для выбора соответствующего количества компонентов для гауссовой модели смеси. Например,

Гауссовы модели смесей требуют, чтобы перед подгонкой к данным было указано количество компонентов. Для многих приложений может быть трудно узнать соответствующее количество компонентов. В этом примере используется статистика соответствия AIC, которая помогает выбрать наиболее подходящую модель гауссовой смеси для различных компонентов.

Создание данных из смеси двух двумерных гауссовых распределений.

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

Предполагая, что вы не знаете базовых значений параметров, графики рассеяния предполагают:

Поместите двухкомпонентную модель гауссовой смеси. На основе проверки графика рассеяния укажите, что ковариационные матрицы являются диагональными. Распечатайте окончательную итерацию и статистику средств к существованию в командном окне, передав statset структура как значение Options аргумент пары имя-значение.

options = statset('Display','final');
GMModel = fitgmdist(X,2,'CovarianceType','diagonal','Options',options);

11 iterations, log-likelihood = -4787.38

GMModel является встроенным gmdistribution модель.

Осмотрите АПК на различное количество компонентов.

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 возникает, когда программное обеспечение подходит под двухкомпонентную модель гауссовой смеси.

Оценки параметров модели гауссовой смеси могут изменяться с различными исходными значениями. В этом примере показано, как управлять начальными значениями при подгонке гауссовых моделей смесей с помощью fitgmdist.

Загрузите набор данных радужки Фишера. Используйте длины и ширину лепестков в качестве предикторов.

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

Поместите модель смеси по Гауссу в данные, используя начальные значения по умолчанию. Существует три вида радужки, поэтому укажите k = 3 компонента.

rng(10); % For reproducibility
GMModel1 = fitgmdist(X,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

По контурам, 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 кажется, что лучше всего различают виды радужки.