之前學習了貝葉斯分類器的構造和使用，其中核心的部分是**得到事件的先驗概率并計算出后驗概率**?，而事實上在實際使用中，很多時候無法得到這些完整的信息，因此我們需要使用另外一個重要的工具——參數估計。 參數估計是在已知系統模型結構時，用系統的輸入和輸出數據計算系統模型參數的過程。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](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 ~~~