文章轉自：http://www.kancloud.cn/digest/prandmethod/102846

        之前學習了貝葉斯分類器的構造和使用，其中核心的部分是得到事件的先驗概率並計算出後驗概率 ，而事實上在實際使用中，很多時候無法得到這些完整的資訊，因此我們需要使用另外一個重要的工具——引數估計。

引數估計是在已知系統模型結構時，用系統的輸入和輸出資料計算系統模型引數的過程。18世紀末德國數學家C.F.高斯首先提出引數估計的方法，他用最小二乘法計算天體執行的軌道。20世紀60年代，隨著電子計算機的普及，引數估計有了飛速的發展。引數估計有多種方法，有最小二乘法、極大似然法、極大驗後法、最小風險法和極小化極大熵法等。在一定條件下，後面三個方法都與極大似然法相同。最基本的方法是最小二乘法和極大似然法。這裡使用MATLAB實現幾種經典的方法並研究其估計特性。

目前理論部分尚在研究當中，無法解釋得很清楚，這是一個長期的學習過程。

一、最大似然估計法

最大似然估計是一種統計方法，它用來求一個樣本集的相關概率密度函式的引數。這個方法最早是遺傳學家以及統計學家羅納德•費雪爵士在 1912 年至1922 年間開始使用的。這種方法的基本思想是：當從模型總體隨機抽取n組樣本觀測值後，最合理的引數估計量應該使得從模型中抽取該n組樣本觀測值的概率最大，而不是像最小二乘估計法旨在得到使得模型能最好地擬合樣本資料的引數估計量。

最大似然估計的優點主要有： 
1. 隨著樣本數量的增加，收斂性變好； 
2. 比任何其他的迭代技術都簡單，適合實用。

最大四然估計的一般原理：

假設有c個類，且：

在正態分佈時有：

假設資料集合D劃分成的類為

我們可以將資料集合D劃分成互不相交的樣本子集合，每一個子集合中的樣本屬於同一類。對每一個數據集合Di ，單獨估計自己的 ，這樣只需估計其引數即可得到分佈函式。

這樣，假設集合D包含同一類n個樣本, x1, x2,…, xn，且這些樣本是獨立抽樣得到的，則： 

P(D | θ)稱為樣本集合Ｄ的似然函式。引數θ 的最大似然估計是通過使定義的P(D | θ) 最大化得到的θ的值，使得實際觀測到的樣本集合具有最大的概率。在實際中，可以對似然函式進行對數運算，為此定義對數似然函式如下：

這樣有：

其導數定義為：

（注：由於自然對數函式是單調增函式，因此對數似然函式和似然函式的極值點的位置相同！）

這裡我們要尋找最優解，這是最大似然估計的基本思想。最優解的必要條件是：

若θ由p個引數組成，則上式代表p個方程組成的方程組。下一步是求取似然函式的極值點：

此時問題可以重新表述為：求對數似然函式的極值點引數θ，使對數似然函式取得最大值：

根據以上公式，在樣本點服從正態分佈，且均值和協方差矩陣未知的情況下，我們可以得到以下兩個公式，分別計算估計均值和協方差矩陣。

二、貝葉斯估計法

在最大似然估計中，假設引數θ未知，但事實上它是確定的量。在而在貝葉斯估計法中，引數θ視為一個隨機變數。後驗概率P(x| D) 的計算是貝葉斯學習的最終目的，其核心問題是給定樣本集合D，計算P(θ| D)。方法二和方法三就屬於貝葉斯估計法。（公式的推導有待進一步研究）

離散情況：

三、估計方法的公式總結

這裡假設樣本點服從正態分佈，且均值和協方差矩陣未知的情況下，方法一計算估計均值和協方差矩陣：

方法二計算估計均值和協方差矩陣：

方法三計算估計均值和協方差矩陣：

其中，λ 是一個小於1的正數，一般可取λ=0.01。可以看到，方法二和方法三均使用迭代的方法求解每一步的均值向量和協方差矩陣。

四、估計誤差

在方法二、三中，在第k步的估計誤差可以由下面的公式計算：

上面的是向量和矩陣範數，一般可以取歐氏範數的形式。

五、實驗步驟與討論

對不同維數下的高斯概率密度模型，用最大似然估計等方法對其引數進行估計。假設有概率密度為高斯分佈p(x) ~N(μ,∑)的樣本集合S={xk,k = 1,2,…,n}，可以使用以上介紹的三種方法進行均值向量和協方差矩陣的引數估計，實驗步驟主要包括以下幾個部分：

1)參考:http://blog.csdn.net/liyuefeilong/article/details/44247589 中的程式，呼叫模式類生成函式MixGaussian，要求給定一組均值（均值向量）和方差（協方差矩陣），生成一維、二維或三維的高斯分佈的散點，每一組資料集合的樣本個數為100。

2)利用上面介紹的方法一、方法二和方法三，分別估計這三組資料的均值（均值向量）和方差（協方差矩陣）。

3)利用上面介紹的誤差計算公式，計算方法一、方法二和方法三的估計誤差。對於方法二和方法三，計算每一次迭代估計的誤差值，並繪製出誤差曲線圖。

這裡生成三維散點圖，以下是三種方法的估計結果與誤差曲線圖：

方法一：

方法二：

方法三：

結論一： 觀察誤差曲線圖，方法三在估計樣本點個數少的時候比方法二的均值誤差更小，但方法二的均值誤差曲線總體呈單調遞減的走勢，而方法三正好相反。在觀察誤差曲線時，無法找出更多可靠的規律。考慮到這裡樣本點只有100個，因此增加樣本點個數進行實驗。

4)在(1)中，將每一個數據集合的樣本個數變為10000，重新生成資料集合，重複(2) (3)中的實驗。

方法一：

方法二：

方法三：

結論二：可以看到，使用方法二進行迭代求取均值和協方差矩陣的方法的誤差曲線隨著樣本點數的增多而區域平緩；相比之下，方法三的誤差曲線變化幅度較大，而且在樣本個數較多的情況下誤差要大於方法二。

5)該部分是論證方法三中λ的取值對估計誤差的影響。這裡選取不同的λ值，重複(2) (3) (4)中的實驗，觀察結果的異同，並加以解釋分析。

100資料lamda比較：

10000資料lamda比較：

結論三：結合輸出結果可以發現，方法三中λ的取值對均值向量估計和協方差矩陣估計的誤差有很大影響，過大的λ取值會導致更大的估計誤差，因此通常λ取0.01。

總結：由於實驗採用了有限樣本估計帶來了一定的估計誤差，這一誤差的影響可以通過增加樣本個數的方法來減小。在理論上，貝葉斯估計方法有很強的理論和演算法基礎，但實際中，最大似然估計演算法更加簡便，而且設計出的分類器效能幾乎與貝葉斯方法得到的結果相差無幾。若碰到樣本點個數很少時，可以嘗試使用方法三進行估計。

最後給出幾個主要函式的定義，只需簡單設定引數就可以呼叫這些函式：

% 方法一進行引數估計，輸出估計均值和協方差矩陣
function [mu,sigma] = ParameterEstimation1(r)

[x,y] = size(r);
% 散點的均值計算
mu = sum(r)/x;
% 散點的協方差矩陣計算
a = 0;
for i = 1:x
    a = a + (r(i,:,:) - mu)' * (r(i,:,:) - mu);
end
sigma = a/x;

disp(['使用第一種方法估計的',num2str(y),'維均值向量為：']);
disp(mu);
disp(['使用第一種方法估計的',num2str(y),'維協方差矩陣為：']);
disp(sigma);

if y == 1
    k=-3:0.1:x;
    x=mu + 0.*k; % 產生一個和y長度相同的x陣列
    plot(x,k,'g-'); 
else if y == 2
    figure;
    plot(r(:,1),r(:,2),'r.');
    hold on;
    plot(mu(1,1),mu(1,2),'bo');
else if y == 3
    figure;
    plot3(r(:,1),r(:,2),r(:,3),'r.');
    hold on;
    plot3(mu(1,1),mu(1,2),mu(1,3),'bo');
    end
end
end
title('使用方法一得出散點圖與估算均值（藍色圈圈）');

% 方法二進行引數估計，輸出估計均值和協方差矩陣
function [mu,sigma] = ParameterEstimation2(r)

[x,y] = size(r);

% 套用公式前均值和協方差矩陣的初值
mu = zeros(1,y);
sigma = zeros(y,y);

for i = 1:x
mu(:,:,i+1) = mu(:,:,i)*(i-1)/i + r(i,:)/i;
sigma(:,:,i+1) = sigma(:,:,i)*(i-1)/i + (r(i,:) - mu(:,:,i+1))' * (r(i,:) - mu(:,:,i+1)) ./i;
end
disp(['使用第二種方法估計的',num2str(y),'維均值向量為：']);
disp(mu(1,:,x+1));
disp(['使用第二種方法估計的',num2str(y),'維協方差矩陣為：']);
disp(sigma(:,:,x+1));

if y == 1
    k=-3:0.1:x;
    x=mu(1,1,x+1) + 0.*k; % 產生一個和y長度相同的x陣列
    plot(x,k,'g-'); 
else if y == 2
    figure;
    plot(r(:,1),r(:,2),'r.');
    hold on;
    plot(mu(1,1,x+1),mu(1,2,x+1),'bo');
else if y == 3
    figure;
    plot3(r(:,1),r(:,2),r(:,3),'r.');
    hold on;
    plot3(mu(1,1,x+1),mu(1,2,x+1),mu(1,3,x+1),'bo');
    end
end
end
title('使用方法二得出散點圖與估算均值（藍色圈圈）');

% 方法三進行引數估計，輸出估計均值和協方差矩陣
% 輸入不同維度的高斯概率密度模型，估計均值和協方差矩陣
% lamda為方法三中的調節引數，預設為0.01
function [mu,sigma] = ParameterEstimation3(r, lamda)

% lamda的預設值為0.01
default_lamda = 0.01;   
% 如果只輸入r而沒有設定lamda，則lamda選用預設值
if nargin == 1
    lamda = default_lamda;
end
[x,y] = size(r);

% 套用公式前均值和協方差矩陣的初值
mu = zeros(1,y);
sigma = zeros(y,y);

for i = 1:x
mu(:,:,i+1) = mu(:,:,i)*(1 - lamda) + r(i,:) * lamda;
sigma(:,:,i+1) = sigma(:,:,i)*(1 - lamda) + (r(i,:) - mu(:,:,i+1))' * (r(i,:) - mu(:,:,i+1)) .*lamda;
end
disp(['使用第三種方法估計的',num2str(y),'維均值向量為：']);
disp(mu(1,:,x+1));
disp(['使用第三種方法估計的',num2str(y),'維協方差矩陣為：']);
disp(sigma(:,:,x+1));

if y == 1
    k=-3:0.1:x;
    x=mu(1,1,x+1) + 0.*k; % 產生一個和y長度相同的x陣列
    plot(x,k,'g-'); 
else if y == 2
    figure;
    plot(r(:,1),r(:,2),'r.');
    hold on;
    plot(mu(1,1,x+1),mu(1,2,x+1),'bo');
else if y == 3
    figure;    
    plot3(r(:,1),r(:,2),r(:,3),'r.');
    hold on;
    plot3(mu(1,1,x+1),mu(1,2,x+1),mu(1,3,x+1),'bo');
    end
    end
end
title('使用方法三得出散點圖與估算均值（藍色圈圈）');

% 誤差估計函式
function ErrorEstimate(mu, sigma, Est_mu, Est_sigma)

[x,y,z] = size(Est_mu);
if z == 1
    Emu = norm((Est_mu - mu),2);
    Esigma = norm((Est_sigma - sigma),2);
    fprintf(['使用方法一估計的',num2str(y),'維高斯分佈的均值誤差為：',num2str(Emu),'\n']);
    fprintf(['\n使用方法一估計的',num2str(y),'維高斯分佈的協方差誤差為：',num2str(Esigma),'\n\n']);
else  
for i = 1:z-1
    Emu(i) = norm((Est_mu(:,:,i+1) - mu),2);
    Esigma(i) = norm((Est_sigma(:,:,i+1) - sigma),2);
end
x = 1:1:z-1;
figure;
subplot(1,2,1);
plot(x,Emu,'r.-');
title('均值向量估計誤差曲線圖');
grid on;
subplot(1,2,2),plot(x,Esigma,'r.-');
title('協方差矩陣估計誤差曲線圖');
grid on;
end