下面用MATLAB實現正態分佈的ML估計

% 二維正態分佈的兩分類問題  （ML估計）

clc;
clear;

% 兩個類別資料的均值向量
Mu = [0 0; 3 3]';
% 協方差矩陣
S1 = 0.8 * eye(2);
S(:, :, 1) = S1;
S(:, :, 2) = S1;
% 先驗概率（類別分佈）
P = [1/3 2/3]';
% 樣本資料規模
% 收斂性：無偏或者漸進無偏，當樣本數目增加時，收斂性質會更好
N = 500;


% 1.生成訓練和測試資料
%{
    生成訓練樣本
    N = 500,  c = 2, d = 2
    μ1=[0, 0]'   μ2=[3, 3]'
    S1=S2=[0.8, 0; 0.8, 0]
    p(w1)=1/3   p(w2)=2/3
%}
randn('seed', 0);
[X_train, Y_train] = generate_gauss_classes(Mu, S, P, N);

figure();
hold on;
class1_data = X_train(:, Y_train==1);
class2_data = X_train(:, Y_train==2);
plot(class1_data(1, :), class1_data(2, :), 'r.');
plot(class2_data(1, :), class2_data(2, :), 'g.');
grid on;
title('訓練樣本');
xlabel('N=500');

%{
    用同樣的方法生成測試樣本
    N = 500,  c = 2, d = 2
    μ1=[0, 0]'   μ2=[3, 3]'
    S1=S2=[0.8, 0; 0.8, 0]
    p(w1)=1/3   p(w2)=2/3
%}
randn('seed', 100);
[X_test, Y_test] = generate_gauss_classes(Mu, S, P, N);
figure();
hold on;
test1_data = X_test(:, Y_test==1);
test2_data = X_test(:, Y_test==2);
plot(test1_data(1, :), test1_data(2, :), 'r.');
plot(test2_data(1, :), test2_data(2, :), 'g.');
grid on;
title('測試樣本');
xlabel('N=500');


% 2.用訓練樣本採用ML方法估計引數
% 各類樣本只包含本類分佈的資訊，也就是說不同類別的引數在函式上是獨立的
[mu1_hat, s1_hat] = gaussian_ML_estimate(class1_data);
[mu2_hat, s2_hat] = gaussian_ML_estimate(class2_data);
mu_hat = [mu1_hat, mu2_hat];
s_hat = (1/2) * (s1_hat + s2_hat);


% 3.用測試樣本和估計出的引數進行分類
% 使用歐式距離進行分類
z_euclidean = euclidean_classifier(mu_hat, X_test);
% 使用貝葉斯方法進行分類
z_bayesian = bayes_classifier(Mu, S, P, X_test);


% 4.計算不同方法分類的誤差
err_euclidean = ( 1-length(find(Y_test == z_euclidean')) /  length(Y_test) );
err_bayesian = ( 1-length(find(Y_test == z_bayesian')) /  length(Y_test) );
% 輸出資訊
disp(['基於歐式距離分類的誤分率：', num2str(err_euclidean)]);
disp(['基於最小錯誤率貝葉斯分類的誤分率：', num2str(err_bayesian)]);


% 畫圖展示
figure();
hold on;
z_euclidean = transpose(z_euclidean);
o = 1;
q = 1;
for i = 1:size(X_test, 2)
    if Y_test(i) ~= z_euclidean(i)
        plot(X_test(1,i), X_test(2,i), 'bo');
    elseif z_euclidean(i)==1
        euclidean_classifier_results1(:, o) = X_test(:, i);
        o = o+1;
    elseif z_euclidean(i)==2
        euclidean_classifier_results2(:, q) = X_test(:, i);
        q = q+1;
    end
end
plot(euclidean_classifier_results1(1, :), euclidean_classifier_results1(2, :), 'r.');
plot(euclidean_classifier_results2(1, :), euclidean_classifier_results2(2, :), 'g.');
title(['基於歐式距離分類，誤分率為：', num2str(err_euclidean)]);
grid on;

figure();
hold on;
z_bayesian = transpose(z_bayesian);
o = 1;
q = 1;
for i = 1:size(X_test, 2)
    if Y_test(i) ~= z_bayesian(i)
        plot(X_test(1,i), X_test(2,i), 'bo');
    elseif z_bayesian(i)==1
        bayesian_classifier_results1(:, o) = X_test(:, i);
        o = o+1;
    elseif z_bayesian(i)==2
        bayesian_classifier_results2(:, q) = X_test(:, i);
        q = q+1;
    end
end
plot(bayesian_classifier_results1(1, :), bayesian_classifier_results1(2, :), 'r.');
plot(bayesian_classifier_results2(1, :), bayesian_classifier_results2(2, :), 'g.');
title(['基於最小錯誤率的貝葉斯決策分類，誤分率為：', num2str(err_bayesian)]);
grid on;
 

function [ data, C ] = generate_gauss_classes( M, S, P, N )
%{
    函式功能：
        生成樣本資料，符合正態分佈

    引數說明：
        M：資料的均值向量
        S：資料的協方差矩陣
        P：各類樣本的先驗概率，即類別分佈
        N：樣本規模

    函式返回
        data：樣本資料（2*N維矩陣）
        C：樣本資料的類別資訊
%}

[~, c] = size(M);
data = [];
C = [];

for j = 1:c
    % z = mvnrnd(mu,sigma,n);
    % 產生多維正態隨機數，mu為期望向量，sigma為協方差矩陣，n為規模。
    % fix 函式向零方向取整
    t = mvnrnd(M(:,j), S(:,:,j), fix(P(j)*N))';
    
    data = [data t];
    C = [C ones(1, fix(P(j) * N)) * j];
end

end
 

function [ m_hat, s_hat ] = gaussian_ML_estimate( X )
%{
    函式功能：
        樣本正態分佈的最大似然估計

    引數說明：
        X：訓練樣本

    函式返回：
        m_hat：樣本由極大似然估計得出的正態分佈引數，均值
        s_hat：樣本由極大似然估計得出的正態分佈引數，方差
%}

% 樣本規模
[~, N] = size(X);
% 正態分佈樣本總體的未知均值μ的極大似然估計就是訓練樣本的算術平均
m_hat = (1/N) * sum(transpose(X))';

% 正態分佈中的協方差陣Σ的最大似然估計量等於N個矩陣的算術平均值
s_hat = zeros(1);
for k = 1:N
    s_hat = s_hat + (X(:, k)-m_hat) * (X(:, k)-m_hat)';
end
s_hat = (1/N)*s_hat;
end

% 詳細的計算過程推導可以參考前一篇部落格：極大似然估計詳解。
 

基於歐式距離的分類：

function [ z ] = euclidean_classifier( m, X )
%{
    函式功能：
        利用歐式距離對測試資料進行分類

    引數說明：
        m：資料的均值，由ML對訓練資料，引數估計得到
        X：我們需要測試的資料

    函式返回：
        z：資料所屬的分類
%}

[~, c] = size(m);
[~, n] = size(X);

z = zeros(n, 1);
de = zeros(c, 1);
for i = 1:n
    for j = 1:c
        de(j) = sqrt( (X(:,i)-m(:,j))' * (X(:,i)-m(:,j)) );
    end
    [~, z(i)] = min(de);
end

end

function [ z ] = bayes_classifier( m, S, P, X )
%{
    函式功能：
        利用基於最小錯誤率的貝葉斯對測試資料進行分類

    引數說明：
        m：資料的均值
        S：資料的協方差
        P：資料類別分佈概率
        X：我們需要測試的資料

    函式返回：
        z：資料所屬的分類
%}

[~, c] = size(m);
[~, n] = size(X);

z = zeros(n, 1);
t = zeros(c, 1);
for i = 1:n
    for j = 1:c
        t(j) = P(j) * comp_gauss_dens_val( m(:,j), S(:,:,j), X(:,i) );
    end
    [~, z(i)] = max(t);
end

end

function [ z ] = comp_gauss_dens_val( m, s, x )
%{
    函式功能：
        計算高斯分佈N（m, s），在某一個特定點的值

    引數說明：
        m：資料的均值
        s：資料的協方差
        x：我們需要計算的資料點

    函式返回：
        z：高斯分佈在x出的值
%}

z = ( 1/( (2*pi)^(1/2)*det(s)^0.5 ) ) * exp( -0.5*(x-m)'*inv(s)*(x-m) );

end