深度學習之最大似然估計

參考鏈接

參考鏈接1

參考連接2

參考鏈接3

一、定義

 

 

 二、知識解讀

　　極大似然估計，通俗理解來說，就是利用已知的樣本結果信息，反推最具有可能（最大概率）導致這些樣本結果出現的模型參數值！

　　換句話說，極大似然估計提供了一種給定觀察數據來評估模型參數的方法，即：“模型已定，參數未知”。

　　可能有小伙伴就要說了，還是有點抽象呀。我們這樣想，一當模型滿足某個分布，它的參數值我通過極大似然估計法求出來的話。比如正態分布中公式如下：

 

　　如果我通過極大似然估計，得到模型中參數和的值，那么這個模型的均值和方差以及其它所有的信息我們是不是就知道了呢。確實是這樣的。極大似然估計中采樣需滿足一個重要的假設，就是所有的采樣都是獨立同分布的。

　　下面我通過倆個例子來幫助理解一下最大似然估計

　　但是首先看一下似然函數  的理解：

 

　　對於這個函數：  輸入有兩個：x表示某一個具體的數據；  表示模型的參數

　　如果  是已知確定的，  是變量，這個函數叫做概率函數(probability function)，它描述對於不同的樣本點  ，其出現概率是多少。

　　如果  是已知確定的，  是變量，這個函數叫做似然函數(likelihood function), 它描述對於不同的模型參數，出現  這個樣本點的概率是多少。

　　這有點像“一菜兩吃”的意思。其實這樣的形式我們以前也不是沒遇到過。例如，  , 即x的y次方。如果x是已知確定的(例如x=2)，這就是  , 這是指數函數。 如果y是已知確定的(例如y=2)，這就是  ，這是二次函數。同一個數學形式，從不同的變量角度觀察，可以有不同的名字。

　　這么說應該清楚了吧？ 如果還沒講清楚，別急，下文會有具體例子。

　　現在真要先講講MLE了。。

　　例子一

　　別人博客的一個例子。

　　假如有一個罐子，里面有黑白兩種顏色的球，數目多少不知，兩種顏色的比例也不知。我 們想知道罐中白球和黑球的比例，但我們不能把罐中的球全部拿出來數。現在我們可以每次任意從已經搖勻的罐中拿一個球出來，記錄球的顏色，然后把拿出來的球 再放回罐中。這個過程可以重復，我們可以用記錄的球的顏色來估計罐中黑白球的比例。假如在前面的一百次重復記錄中，有七十次是白球，請問罐中白球所占的比例最有可能是多少？

　　很多人馬上就有答案了：70%。而其后的理論支撐是什么呢？

　　我們假設罐中白球的比例是p，那么黑球的比例就是1-p。因為每抽一個球出來，在記錄顏色之后，我們把抽出的球放回了罐中並搖勻，所以每次抽出來的球的顏 色服從同一獨立分布。

　　這里我們把一次抽出來球的顏色稱為一次抽樣。題目中在一百次抽樣中，七十次是白球的,三十次為黑球事件的概率是P(樣本結果|Model)。

　　如果第一次抽象的結果記為x1,第二次抽樣的結果記為x2....那么樣本結果為(x1,x2.....,x100)。這樣，我們可以得到如下表達式：

　　P(樣本結果|Model)

　　　　= P(x1,x2,…,x100|Model)

　　　　= P(x1|Mel)P(x2|M)…P(x100|M)

　　　　= p^70(1-p)^30.

　　好的，我們已經有了觀察樣本結果出現的概率表達式了。那么我們要求的模型的參數，也就是求的式中的p。

　　那么我們怎么來求這個p呢？

　　不同的p，直接導致P（樣本結果|Model）的不同。

　　好的，我們的p實際上是有無數多種分布的。如下：

 

　　那么求出 p^70(1-p)^30為 7.8 * 10^(-31)

 

　　p的分布也可以是如下：

 

　　那么也可以求出p^70(1-p)^30為2.95* 10^(-27)

　　那么問題來了，既然有無數種分布可以選擇，極大似然估計應該按照什么原則去選取這個分布呢？

　　答：采取的方法是讓這個樣本結果出現的可能性最大，也就是使得p^70(1-p)^30值最大，那么我們就可以看成是p的方程，求導即可！

　　那么既然事情已經發生了，為什么不讓這個出現的結果的可能性最大呢？這也就是最大似然估計的核心。

　　我們想辦法讓觀察樣本出現的概率最大，轉換為數學問題就是使得：

　　p^70(1-p)^30最大，這太簡單了，未知數只有一個p，我們令其導數為0，即可求出p為70%，與我們一開始認為的70%是一致的。其中蘊含着我們的數學思想在里面。

　　使用matlab查看其最大值

clear
clc
x = 0:0.1:1;
y = (x.^70).*((1-x).^30);
plot(x,y)

　　

 

 

 

 

　　例子二

　　假設我們要統計全國人民的年均收入，首先假設這個收入服從服從正態分布，但是該分布的均值與方差未知。我們沒有人力與物力去統計全國每個人的收入。我們國家有10幾億人口呢？那么豈不是沒有辦法了？

　　不不不，有了極大似然估計之后，我們可以采用嘛！我們比如選取一個城市，或者一個鄉鎮的人口收入，作為我們的觀察樣本結果。然后通過最大似然估計來獲取上述假設中的正態分布的參數。

　　有了參數的結果后，我們就可以知道該正態分布的期望和方差了。也就是我們通過了一個小樣本的采樣，反過來知道了全國人民年收入的一系列重要的數學指標量！

　　那么我們就知道了極大似然估計的核心關鍵就是對於一些情況，樣本太多，無法得出分布的參數值，可以采樣小樣本后，利用極大似然估計獲取假設中分布的參數值。

　　例子三

 

 

  三、求解過程

　　由於樣本集中的樣本都是獨立同分布，可以只考慮一類樣本集D，來估計參數向量θ。記已知的樣本集為：

 

 

 

 

        似然函數（linkehood function）：聯合概率密度函數稱為相對於的θ的似然函數。

 

 

 

 

        如果是參數空間中能使似然函數最大的θ值，則應該是“最可能”的參數值，那么就是θ的極大似然估計量。它是樣本集的函數，記作：

 

 

 

求解極大似然函數

        ML估計：求使得出現該組樣本的概率最大的θ值。

 

         實際中為了便於分析，定義了對數似然函數：

  1. 未知參數只有一個（θ為標量）

        在似然函數滿足連續、可微的正則條件下，極大似然估計量是下面微分方程的解：

 

        2.未知參數有多個（θ為向量）

        則θ可表示為具有S個分量的未知向量：

 

         記梯度算子：

 

 

 

 

         若似然函數滿足連續可導的條件，則最大似然估計量就是如下方程的解。

 

 

 

 

         方程的解只是一個估計值，只有在樣本數趨於無限多的時候，它才會接近於真實值。

 

極大似然估計的例子
        例1：設樣本服從正態分布，則似然函數為：

        它的對數：

 

        求導，得方程組：

        聯合解得：

 

 

 

 

        似然方程有唯一解，而且它一定是最大值點，這是因為當時，非負函數。於是U和的極大似然估計為。

 

        例2：設樣本服從均勻分布[a, b]。則X的概率密度函數：

 

        對樣本

 

 

 

        很顯然，L(a,b)作為a和b的二元函數是不連續的，這時不能用導數來求解。而必須從極大似然估計的定義出發，求L(a,b)的最大值，為使L(a,b)達到最大，b-a應該盡可能地小，但b又不能小於，否則，L(a,b)=0。類似地a不能大過，因此，a和b的極大似然估計：

 

 

 

總結
        求最大似然估計量的一般步驟：

        （1）寫出似然函數；

        （2）對似然函數取對數，並整理；

        （3）求導數；

        （4）解似然方程。

        最大似然估計的特點：

        1.比其他估計方法更加簡單；

        2.收斂性：無偏或者漸近無偏，當樣本數目增加時，收斂性質會更好；

        3.如果假設的類條件概率模型正確，則通常能獲得較好的結果。但如果假設模型出現偏差，將導致非常差的估計結果。
Matlab實例

　　實例一

 

 

 

x=[1.2 3.5 4.2 0.8 1.4 3.1 4.8 0.9];
u=mle(x)

　　實例二

　　下面用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

　　正態分布的ML估計（對訓練樣本）：

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 ] = 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