1. 程式人生 > >機器學習演算法(1)——極大似然估計與EM演算法

機器學習演算法(1)——極大似然估計與EM演算法

極大似然估計

在講解極大似然估計前,需要先介紹貝葉斯分類:

貝葉斯決策:

      首先來看貝葉斯分類,經典的貝葉斯公式:

                                                               

             其中:p(w)為先驗概率,表示每種類別分佈的概率;P(x|w)是條件概率,表示在某種類別前提下,某件事發生的概率;而P(w|x)為後驗概率,表示某事發生了,並且它屬於某一類別的概率,有了這個後驗概率,我們就可以對樣本進行分類。後驗概率越大,說明某事物屬於這個類別的可能性越大,我們越有理由把它歸到這個類別下。

        我們來看一個直觀的例子:已知:某個班級裡男生中學文科的概率為1/2,女生中學文科的概率為2/3,並且該班級中男女比例為2:1,問題:若你在班級遇到一個學文科的人,請問他的性別為男或女的概率分別為多少?問題就是,某事發生了,它屬於某一類別的概率是多少?即後驗概率。
設:w_{1}=man,w_{2}=girl,x=liberal

  由已知得:

    

由於男女選擇學文科相互獨立,所以選擇學文科的概率為:

                                P(x)=P(x|w_{1},w_{2})=P(x|w_{1})*P(w_{1})+P(x|w_{2})*P(w_{2})=5/9

考慮性別分類問題,只需要比較後驗概率的大小,取值並不重要,由貝葉斯公式求得:

                                  

     但是在實際問題中,我們能獲得的資料可能只有有限數目的樣本資料,而先驗概率P(w_{i})和條件概率P(x|w_{i})(各類的總體分佈)都是未知的。根據僅有的樣本資料進行分類時,一種可行的辦法是我們需要先對先驗概率和類條件概率進行估計,然後再套用貝葉斯分類器。

        先驗概率P(w_{i})的估計較簡單,1、每個樣本所屬的自然狀態都是已知的(有監督學習);2、依靠經驗;3、用訓練樣本中各類出現的頻率估計。

        類條件概率P(x|w_{i})的估計(非常難),原因包括:概率密度函式包含了一個隨機變數的全部資訊;樣本資料可能不多;特徵向量x的維度可能很大等等。總之要直接估計類條件概率的密度函式很難。解決的辦法就是,把估計完全未知的概率密度轉化為估計引數。這裡就將概率密度估計問題轉化為引數估計問題,極大似然估計就是一種引數估計方法。當然了,概率密度函式的選取很重要,模型正確,在樣本區域無窮時,我們會得到較準確的估計值,如果模型都錯了,那估計半天的引數,肯定也沒啥意義了.

重要前提:訓練樣本的分佈能代表樣本的真實分佈。每個樣本集中的樣本都是所謂獨立同分布的隨機變數 (iid條件),且有充分的訓練樣本.

極大似然估計

       舉個別人部落格中的例子,假如有一個罐子,裡面有黑白兩種顏色的球,數目多少不知,兩種顏色的比例也不知。我 們想知道罐中白球和黑球的比例,但我們不能把罐中的球全部拿出來數。現在我們可以每次任意從已經搖勻的罐中拿一個球出來,記錄球的顏色,然後把拿出來的球 再放回罐中。這個過程可以重複,我們可以用記錄的球的顏色來估計罐中黑白球的比例。假如在前面的一百次重複記錄中,有七十次是白球,請問罐中白球所佔的比例最有可能是多少?很多人馬上就有答案了:70%。而其後的理論支撐是什麼呢?

      我們假設罐中白球的比例是p,那麼黑球的比例就是1-p。因為每抽一個球出來,在記錄顏色之後,我們把抽出的球放回了罐中並搖勻,所以每次抽出來的球的顏 色服從同一獨立分佈。這裡我們把一次抽出來球的顏色稱為一次抽樣。題目中在一百次抽樣中,七十次是白球的概率是P(Data | M),這裡Data是所有的資料,M是所給出的模型,表示每次抽出來的球是白色的概率為p。如果第一抽樣的結果記為x1,第二抽樣的結果記為x2... 那麼Data = (x1,x2,…,x100).

  似然函式:P(Data | M)= P(x1,x2,…,x100|M)= P(x1|M)P(x2|M)…P(x100|M)= p^70(1-p)^30;

   p在取什麼值的時候,P(Data |M)的值最大呢?將p^70(1-p)^30對p求導,並其等於零。70p^69(1-p)^30-p^70*30(1-p)^29=0。

解方程可以得到p=0.7。當p=0.7時,P(Data|M)的值最大。這和我們常識中按抽樣中的比例來計算的結果是一樣的。

由上可知最大似然估計的一般求解過程:

  (1) 寫出似然函式;

  (2) 對似然函式取對數,並整理;

  (3) 求導數 ;

  (4) 解似然方程

用MATLAB實現正態分佈的ML估計:from:https://blog.csdn.net/zengxiantao1994/article/details/72793608

EM演算法

好吧我也不知道為什麼將EM演算法,要扯那麼多似然估計。。。。

哈哈當然有用了,還是用例子來說明

1、舉例說明:學生身高問題

  我們需要調查我們學校的男生和女生的身高分佈。 假設你在校園裡隨便找了100個男生和100個女生。他們共200個人。將他們按照性別劃分為兩組,然後先統計抽樣得到的100個男生的身高。假設他們的身高是服從高斯分佈的。但是這個分佈的均值u和方差∂我們不知道,這兩個引數就是我們要估計的。記作θ=[u, ∂]T。

         

問題:已知兩個:(1)樣本服從的分佈模型(2)隨機抽取的樣本 ,需要通過極大似然估計求出的包括:模型的引數。

問題數學化: (1)樣本集X={x1,x2,…,xN} ,N=100 (2)概率密度:p(xi|θ)抽到男生i(的身高)的概率 100個樣本之間獨立同分布,所以我同時抽到這100個男生的概率就是他們各自概率的乘積。就是從分佈是p(x|θ)的總體樣本中抽取到這100個樣本的概率,也就是樣本集X中各個樣本的聯合概率,用下式表示:

這個概率反映了,在概率密度函式的引數是θ時,得到X這組樣本的概率。 需要找到一個引數θ,其對應的似然函式L(θ)最大,也就是說抽到這100個男生(的身高)概率最大。這個叫做θ的最大似然估計量。

然後對似然函式取對數,並整理:

        

然後,求導數,令導數為0,得到似然方程;最後,解似然方程,得到的引數即為所求。

2、傳統EM演算法詳述

       問題:我們抽取的100個男生和100個女生樣本的身高,但是我們不知道抽取的那200個人裡面的每一個人到底是從男生的那個身高分佈裡面抽取的,還是女生的那個身高分佈抽取的(例如樣本x4=1.8m,那麼這個一米八的身高樣本到底是從男生中取出的還是女生中取出的?)。 用數學的語言就是,抽取得到的每個樣本(身高)都不知道是從哪個分佈(是男是女)抽取的。 這個時候,對於每一個樣本,就有兩個東西需要猜測或者估計: (1)這個人是男的還是女的?(2)男生和女生對應的身高的高斯分佈的引數是多少?

EM演算法要解決的問題是: (1)求出每一個樣本屬於哪個分佈 (2)求出每一個分佈對應的引數

    

身高問題使用EM演算法求解步驟:

                           

  步驟:

(1)初始化引數:先初始化男生和女生身高的正態分佈的引數:如均值E_{man1}=1.7,方差D_{man1}=0.1,E_{girl1}=1.6,D_{girl1}=0.2(即初始了兩個概率密度函式)

(2)計算每一個人更可能屬於男生分佈或者女生分佈,哪個概率值大,就認為這個樣本屬於哪一類;(將每個身高分別帶入兩個概率密度函式計算概率,取結果較大值對應的分佈)

(3)通過分為男生的n個人來重新估計男生身高分佈的引數E_{man2}D_{man2}(最大似然估計),女生分佈也按照相同的方式估計出來E_{girl2},D_{girl2},更新分佈。

(4)這時候兩個分佈的概率也變了,然後重複步驟(1)至(3),直到引數不發生變化為止。

4、演算法推導

已知:樣本集X={x(1),…,x(m))},包含m個獨立的樣本;

未知:每個樣本i對應的類別z(i)是未知的(相當於聚類);

輸出:我們需要估計概率模型p(x,z)的引數θ;

目標:找到適合的θ和z讓L(θ)最大

演算法流程:

1)初始化分佈引數θ; 重複以下步驟直到收斂:        

E步驟:根據引數初始值或上一次迭代的模型引數來計算出隱性變數的後驗概率,其實就是隱性變數的期望。作為隱藏變數的現估計值:

 

M步驟:將似然函式最大化以獲得新的引數值:

    

from:https://www.cnblogs.com/Gabby/p/5344658.html

還有一個有趣的例子:from:http://www.cnblogs.com/bigmoyan/p/4550375.html

OpenCV的實現:

EM演算法是一種非監督的學習演算法,它的輸入資料事先不需要進行標註。相反,該演算法從給定的樣本集中,能計算出高斯混和引數的最大似然估計。也能得到每個樣本對應的標註值,類似於kmeans聚類(輸入樣本資料,輸出樣本資料的標註)。實際上,高斯混和模型GMM和kmeans都是EM演算法的應用。

在opencv3.0中,EM演算法的函式是trainEM,函式原型為:

//EM步驟
bool trainEM(InputArray samples, OutputArray logLikelihoods=noArray(),OutputArray labels=noArray(),OutputArray probs=noArray())

//E步驟
bool trainE(InputArray samples, InputArray means0,
                        InputArray covs0=noArray(),
                        InputArray weights0=noArray(),
                        OutputArray logLikelihoods=noArray(),
                        OutputArray labels=noArray(),
                        OutputArray probs=noArray())
//M步驟
bool trainM(InputArray samples, InputArray probs0,
                        OutputArray logLikelihoods=noArray(),
                        OutputArray labels=noArray(),
                        OutputArray probs=noArray())

輸入四個引數:

 samples: 輸入的樣本,一個單通道的矩陣。從這個樣本中,進行高斯混和模型估計。

logLikelihoods: 可選項,輸出一個矩陣,裡面包含每個樣本的似然對數值。

labels: 可選項,輸出每個樣本對應的標註。

probs: 可選項,輸出一個矩陣,裡面包含每個隱性變數的後驗概率

        這個函式沒有輸入引數的初始化值,是因為它會自動執行kmeans演算法,將kmeans演算法得到的結果作為引數初始化。這個trainEM函式實際把E步驟和M步驟都包含進去了,我們也可以對兩個步驟分開執行,OPENCV3.0中也提供了分別執行的函式trainE、trainM。

         trainEM函式的功能和kmeans差不多,都是實現自動聚類,輸出每個樣本對應的標註值。但它比kmeans還多出一個功能,就是它還能起到訓練分類器的作用,用於後續新樣本的預測。預測函式原型為:

Vec2d predict2(InputArray sample, OutputArray probs) const

sample: 待測樣本

probs : 和上面一樣,一個可選的輸出值,包含每個隱性變數的後驗概率。

返回一個Vec2d型別的向量,包括兩個元素,第一個元素為樣本的似然對數值,第二個元素為最大可能混和分量的索引值

聚類加預測:

#include "stdafx.h"
#include "opencv2/opencv.hpp"
#include <iostream>
using namespace std;
using namespace cv;
using namespace cv::ml;

//使用EM演算法實現樣本的聚類及預測
int main()
{
    const int N = 4;    //分成4類
    const int N1 = (int)sqrt((double)N);
    //定義四種顏色,每一類用一種顏色表示
    const Scalar colors[] =
    {
        Scalar(0, 0, 255), Scalar(0, 255, 0),
        Scalar(0, 255, 255), Scalar(255, 255, 0)
    };

    int i, j;
    int nsamples = 100;   //100個樣本點
    Mat samples(nsamples, 2, CV_32FC1);  //樣本矩陣,100行2列,即100個座標點    
    Mat img = Mat::zeros(Size(500, 500), CV_8UC3);  //待測資料,每一個座標點為一個待測資料
    samples = samples.reshape(2, 0);

    //迴圈生成四個類別樣本資料,共樣本100個,每類樣本25個
    for (i = 0; i < N; i++)
    {
        
        Mat samples_part = samples.rowRange(i*nsamples / N, (i + 1)*nsamples / N);

        //設定均值
        Scalar mean(((i%N1) + 1)*img.rows / (N1 + 1),
            ((i / N1) + 1)*img.rows / (N1 + 1));
        //設定標準差
        Scalar sigma(30, 30);
        randn(samples_part, mean, sigma);  //根據均值和標準差,隨機生成25個正態分佈座標點作為樣本
    }
    samples = samples.reshape(1, 0);
    // 訓練分類器
    Mat labels;  //標註,不需要事先知道
    Ptr<EM> em_model = EM::create();
    em_model->setClustersNumber(N);
    em_model->setCovarianceMatrixType(EM::COV_MAT_SPHERICAL);
    em_model->setTermCriteria(TermCriteria(TermCriteria::COUNT + TermCriteria::EPS, 300, 0.1));
    em_model->trainEM(samples, noArray(), labels, noArray());

    //對每個座標點進行分類,並根據類別用不同的顏色畫出
    Mat sample(1, 2, CV_32FC1);
    for (i = 0; i < img.rows; i++)
    {
        for (j = 0; j < img.cols; j++)
        {
            sample.at<float>(0) = (float)j;
            sample.at<float>(1) = (float)i;
            //predict2返回的是double值,用cvRound進行四捨五入得到整型
            //此處返回的是兩個值Vec2d,取第二個值作為樣本標註
            int response = cvRound(em_model->predict2(sample, noArray())[1]);
            Scalar c = colors[response];  //為不同類別設定顏色
            circle(img, Point(j, i), 1, c*0.75, FILLED);
        }
    }

    //畫出樣本點
    for (i = 0; i < nsamples; i++)
    {
        Point pt(cvRound(samples.at<float>(i, 0)), cvRound(samples.at<float>(i, 1)));
        circle(img, pt, 2, colors[labels.at<int>(i)], FILLED);
    }

    imshow("EM聚類結果", img);
    waitKey(0);

    return 0;
}

結果:

只聚類無預測:

Vec3b colorTab[] =
    {
        Vec3b(0, 0, 255),
        Vec3b(0, 255, 0),
        Vec3b(255, 100, 100),
        Vec3b(255, 0, 255),
        Vec3b(0, 255, 255)
    };//建立顏色索引


int N =3;  //聚成3類
    Ptr<EM> em_model = EM::create();
    em_model->setClustersNumber(N);
    em_model->setCovarianceMatrixType(EM::COV_MAT_SPHERICAL);
    em_model->setTermCriteria(TermCriteria(TermCriteria::COUNT + TermCriteria::EPS, 300, 0.1));
    em_model->trainEM(data, noArray(), labels, noArray());

    int n = 0;
    //顯示聚類結果,不同的類別用不同的顏色顯示
    for (int i = 0; i < pic.rows; i++)
    for (int j = 0; j < pic.cols; j++)
    {
        int clusterIdx = labels.at<int>(n);//label存放了該樣本屬於哪個型別的索引
        pic.at<Vec3b>(i, j) = colorTab[clusterIdx];
        n++;
    }
    imshow("pic", pic);

from:https://www.cnblogs.com/denny402/p/5036288.html