捕鱼机遥控器是不是真的管用

書名： Visual C++數字圖像模式識別典型案例詳解
作者名：馮偉興梁洪王臣業編著
本章字數： 3406字
更新時間： 2018-12-31 19:38:57

2.2.2 隱馬爾可夫模型

隱馬爾可夫模型（Hidden Markov Models，HMM）是序列數據處理和統計學習的一種重要概率模型，具有建模簡單、數據計算量小、運行速度快、識別率高等特點，近幾年來已經被成功應用到許多工程任務中。

1.隱馬爾可夫概念

隱馬爾可夫模型（HMM）由兩個序列組成，一個是不可觀察的（隱藏的）狀態變化序列，一個是由該隱藏的狀態所產生的可觀察的符號序列。隱馬爾可夫模型由其初始分布、轉移概率矩陣和發射概率矩陣完全確定。一個建好的隱馬爾可夫模型的運作過程可以這樣描述：由初始分布隨機地產生一個狀態，然后按轉移概率矩陣隨機地從一個狀態轉移到另一個狀態，在每一個狀態按發射概率矩陣隨機地生成一個字符，結果是產生一個狀態序列和一個相應的字符序列。

一般情況下，不僅模型能夠產生的狀態序列不是唯一的，每一個狀態序列能夠產生的字符序列也不是唯一的。模型以一定的概率產生特定的狀態序列，該狀態序列又以一定的概率產生特定的字符序列。在實際問題中，通常字符序列是可觀測的，狀態序列是不可觀測的，是“隱藏的”。

HMM可以表示為λ=(A,S,π)，描述如下：

（a）隱狀態集S。S={S1,S2,S3,…,SN}，N表示狀態數。

（b）觀察符號集V。V ={V1,V2,V3,…,VM}，M表示符號數。

（c）狀態轉換概率矩陣A。A={ai}j

式中

式中，Q=(q1,q2,…,qL)表示狀態序列，其中qi∈S，1≤i≤L，L表示序列的長度。t為當前時刻。aij為從狀態si轉換到s j的概率。

（d）符號釋放概率矩陣E。E={e jk}，其中

式中，O=(o1,o2,…,oL)表示符號序列，其中o j∈V，1≤ j≤L。e jk表示在狀態S j下產生符號Vk的概率。

（e）初始狀態分布π。π={πi}，其中

表示模型起始于狀態Si（即時間為0時的狀態）的概率。

基于HMM的模式識別方法的基本思想：對于一個可用HMM描述的模式識別問題，其可能出現的所有不同的隱馬爾可夫模型λ1,λ2,λ3,…,λn構成一個集合Λ，其中任意一個元素λi為特定的隱馬爾可夫模型，其中包括此特定隱馬爾可夫模型的轉移矩陣、初分布及觀測概率。觀測概率的定義如下：

給定一組樣本觀測值(o1,o2,…,on)，選擇λ?∈Λ，使得

則λ?是觀測(o1,o2,…,on)的最可能模式。

2.隱馬爾可夫模型基本算法

HMM需解決的基本問題：

·概率問題。給定一個HMM（λ=(A,S,π)）和一個觀測符號序列O，求由該模型得到序列O的概率P(O|λ)。

·解碼問題：給定HMM（λ=(A,S,π)）和一個觀測序列O，解碼問題就是尋找最可能的隱狀態序列，即如何選擇產生該觀察序列O的最佳狀態序列Q*，也就是如何最好地解釋符號序列O。

·學習問題。給定一組觀測序列樣本集，如何確定其對應的HMM（λ=(A,S,π)），使該模型能夠最好地描述觀測樣本集，也就是該模型的參數能夠使P(O|λ)最大。

（1）前向后向算法

對于可能性問題，引入前向變量和后向變量進行動態規劃求解。

1）前向算法

前向算法的簡單圖解如圖2-8所示。

圖2-8 前向算法示意圖

前向變量αj(t)表示在t時刻狀態為s j，且序列的前t個符號為o1,o2,…,ot的概率，即

根據前向變量的定義，可以得到：

前向變量的迭代公式：

前向算法的流程：

（a）初始化：

ai(0)=πiei(o 1),1≤i≤N

（b）遞歸：

          t=1,2,3,…,L
              j=1,2,3,…,N
                  a j(t)=ej(ot)∑k(ak(t-1)akj)
              next
          next

（c）終止：

P (o|λ)=∑kαk(L)

（d）算法結束。

2）后向算法

后向算法的簡單圖解如圖2-9所示。

圖2-9 后向算法示意圖

后向變量βj(t)表示在t時刻狀態為s j，且序列從t+1到L個符號為ot+1,ot+2,…,ot+L的概率，即

根據后向變量的定義，可以得到：

后向變量的迭代公式：

后向算法的流程如下：

（a）初始化：

βi (L)=1 1≤i≤N

（b）遞歸：

          t=L-1,…,1
              j= 1,2,…,N
                βj(t)=∑kαjk ek(ot+ 1)βk(t+1)
              next
        next

（c）終止：

Po λ( | )=∑k k k(1) k(1)

（d）算法結束。前向算法和后向算法的時間復雜度為o(N 2L)，空間復雜度為o(NL)。

（2）解碼問題

根據給定的HMM(λ=(A,E,π))和符號序列O=(o1,o2,…,oL)，確定一條具有最大概率的狀態序列，即最優路徑Q*：

應用Viterbi算法求解最優路徑。定義變量δj(t)，它表示由狀態序列q1,q 2,…,qt且qt =s j產生出符號序列o1,o2,…,ot的最大概率，即

δj (t)可以通過迭代計算得到：

從初始狀態開始迭代計算δj(t)，一直到序列末端，然后再回溯得到最優路徑。

Viterbi算法流程如下：

（a）初始化：

δi(1)=πiei(o 1) 1≤i≤N

（b）遞歸：

（c）終止：

（d）回溯：

qt *=φt+1(qt+1)，t=L-1,L-2,…,1

（e）算法結束。

Viterbi算法提供了有效的計算方法來分析HMM中的觀察序列，以此得到最可能的隱狀態序列，并利用遞歸來減小計算負擔。

（3）訓練問題當只有符號序列o1,o2,…,on是給定的時，才要通過訓練學習估計模型的各項參數，使觀測的符號序列能夠“最好”地被描述，即λ能夠使P(o1,o2,…,on|λ)最大。

如果符號序列是相互獨立的，則它們的聯合分布概率即為它們單個概率的乘積：

如果路徑已知，則要統計各個特定的初始狀態，轉移概率以及產生概率。記某種初始狀態、轉移和出現的次數為πi、Aij和Ei(Vk)，則初始狀態概率、轉移概率和產生概率的極大似然概率估計值為

如果路徑未知，常采用Baum-Welch算法來實現對HMM的參數估計。

根據前向變量和后向變量的定義，可以得到：

通過對所有位置以及所有序列的求和，可以得到：

再根據式（2-116）～式（2-118）即可以得到模型的各項參數估計。

Baum-Welch算法流程如下：

定義變量γt(i)=P(qt =si|O,λ)，ξt(i,j)=P(qt =si,qt+1=s j|O,λ)

則

且。

故=時刻t=1在狀態Si的期望值=γ1(i)。

此時有

且，1≤i≤N,1≤ j≤N。

可以通過證明得出結論：且P(o|λ)迭代收斂到局部最佳。

3.隱馬爾可夫模型的C語言實現

下面介紹隱馬爾可夫模型的C語言實現代碼。

動態數組申請的代碼如代碼2-2所示。

代碼2-2動態數組申請代碼

        float ＊vector(long nl,long nh)
        {
            float ＊v;
            v=(float ＊)calloc((unsigned) (nh-nl+1),sizeof(float));
            if (!v) nrerror("allocation failure in vector()");
            return v-nl;
        }
        int ＊ivector(long nl,long nh)
        {
            int ＊v;
            v=(int ＊)calloc((unsigned) (nh-nl+1),sizeof(int));
            if (!v) nrerror("allocation failure in ivector()");
            return v-nl;
        }
        unsigned char ＊cvector(long nl,long nh)
        {
            unsigned char ＊v;
            v=(unsigned char ＊)calloc((unsigned) (nh-nl+1),sizeof(unsigned char));
            if (!v) nrerror("allocation failure in cvector()");
            return v-nl;
        }
        unsigned long ＊lvector(long nl,long nh)
        {
            unsigned long ＊v;
            v=(unsigned long ＊)calloc((unsigned) (nh-nl+1),sizeof(unsigned long));
            if (!v) nrerror("allocation failure in lvector()");
            return v-nl;
        }
        double ＊dvector(long nl,long nh)
        {
            double ＊v;
            v=(double ＊)calloc((unsigned) (nh-nl+1),sizeof(double));
            if (!v) nrerror("allocation failure in dvector()");
            return v-nl;
        }
        float ＊＊matrix(long nrl,long nrh,long ncl,long nch)
        {
            int i;
            float ＊＊m;
            m=(float ＊＊) calloc((unsigned) (nrh-nrl+1),sizeof(float＊));
            if (!m) nrerror("allocation failure 1 in matrix()");
            m -= nrl;
            for(i=nrl;i<=nrh;i++) {
                m[i]=(float ＊) calloc((unsigned) (nch-ncl+1),sizeof(float));
                if (!m[i]) nrerror("allocation failure 2 in matrix()");
                m[i] -= ncl;
            }
            return m;
        }
        double ＊＊dmatrix(long nrl,long nrh,long ncl,long nch)
        {
            int i;
            double ＊＊m;
            m=(double ＊＊) calloc((unsigned) (nrh-nrl+1),sizeof(double＊));
            if (!m) nrerror("allocation failure 1 in dmatrix()");
            m -= nrl;
            for(i=nrl;i<=nrh;i++) {
                m[i]=(double ＊) calloc((unsigned) (nch-ncl+1),sizeof(double));
                if (!m[i]) nrerror("allocation failure 2 in dmatrix()");
                m[i] -= ncl;
            }
            return m;
        }
        int ＊＊imatrix(long nrl,long nrh,long ncl,long nch)
        {
            int i,＊＊m;
            m=(int ＊＊)calloc((unsigned) (nrh-nrl+1),sizeof(int＊));
            if (!m) nrerror("allocation failure 1 in imatrix()");
            m -= nrl;
            for(i=nrl;i<=nrh;i++) {
                m[i]=(int ＊)calloc((unsigned) (nch-ncl+1),sizeof(int));
                if (!m[i]) nrerror("allocation failure 2 in imatrix()");
                m[i] -= ncl;
            }
            return m;
        }
        unsigned char ＊＊cmatrix(long nrl,long nrh,long ncl,long nch)
        {
            int i;
            unsigned char ＊＊m;
            m=(unsigned char ＊＊) calloc((unsigned) (nrh-nrl+1),sizeof(unsigned char＊));
            if (!m) nrerror("allocation failure 1 in dmatrix()");
            m -= nrl;
            for(i=nrl;i<=nrh;i++) {
                m[i]=(unsigned char ＊) calloc((unsigned) (nch-ncl+1),sizeof(unsigned char));
                if (!m[i]) nrerror("allocation failure 2 in dmatrix()");
                m[i] -= ncl;
            }
            return m;
        }
        //求一個二維矩陣的子矩陣
        float ＊＊submatrix(float ＊＊a,long oldrl,long oldrh,long oldcl,long oldch,long
            newrl,long newcl)
        {
            int i,j;
            float ＊＊m;
            m=(float ＊＊) calloc((unsigned) (oldrh-oldrl+1),sizeof(float＊));
            if (!m) nrerror("allocation failure in submatrix()");
            m -= newrl;
            for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+oldcl-newcl;
            return m;
        }
        //將一個實數一維數組轉化為二維數組
        float ＊＊convert_matrix(float ＊a,long nrl,long nrh,long ncl,long nch)
        {
            int i,j,nrow,ncol;
            float ＊＊m;
            nrow=nrh-nrl+1;
            ncol=nch-ncl+1;
            m = (float ＊＊) calloc((unsigned) (nrow),sizeof(float＊));
            if (!m) nrerror("allocation failure in convert_matrix()");
            m -= nrl;
            for(i=0,j=nrl;i<=nrow-1;i++,j++) m[j]=a+ncol＊i-ncl;
            return m;
        }

前向算法的代碼如代碼2-3所示。

代碼2-3前向算法代碼

        //前向算法估計參數
        void Forward(HMM ＊phmm,int T,int ＊O,double ＊＊alpha,double ＊pprob)
        {
            int i,j;
            int t;
            double sum;
            //初始化
            for(i=1;i<=phmm->N;i++)
                alpha[1][i]=phmm->pi[i]＊phmm->B[i][O[1]];
            //遞歸
            for(t=1;t<T;t++){
                for(j=1;j<=phmm->N;j++){
                    sum=0.0;
                for(i=1;j<=phmm->N;i++)
                    sum+=alpha[t][i]＊(phmm->A[i][j]);
                alpha[t+1][j]=sum＊(phmm->B[j][O[t+1]]);
                }
            }
            //終止
            ＊pprob=0.0;
            for(i=1;i<=phmm->N;i++)
                ＊pprob+=alpha[T][i];
        }
        void ForwardWithScale(HMM＊phmm,int T,int ＊O,double ＊＊alpha,double ＊scale,double
            ＊pprob)
        // pprob是對數概率
        {
            int i,j;
            int t;
            double sum;
            //初始化
            scale[1]=0.0;
            for(i=1;i<=phmm->N;i++){
                alpha[1][i]=phmm->pi[i]＊(phmm->B[i][O[1]+1]);
                scale[1]+=alpha[1][i];
            }
            for(i=1;i<=phmm->N;i++)
                alpha[1][i]/=scale[1];
            //遞歸
            for(t=1;t<T;t++){
                scale[t+1]=0.0;
                for(j=1;j<=phmm->N;j++){
                    sum=0.0;
                    for(i=1;i<=phmm->N;i++)
                          sum+=alpha[t][i]＊(phmm->A[i][j]);
                    alpha[t+1][j]=sum＊(phmm->B[j][O[t+1]]);
                    scale[t+1]+=alpha[t+1][j];
                }
                for(j=1;j<=phmm->N;j++)
                    alpha[t+1][j]/=scale[t+1];
            }
            //終止
            ＊pprob=0.0;
            for(t=1;t<=T;t++)
                ＊pprob+=log(scale[t]);
        }

后向算法的代碼如代碼2-4所示。

代碼2-4后向算法代碼

        //后向算法估計參數
        void Backward(HMM ＊phmm, int T, int ＊O, double ＊＊beta, double ＊pprob)
        {
            int i, j;
            int t;
            double sum;
            //初始化
            for (i = 1; i <= phmm->N; i++)
                beta[T][i] = 1.0;
            //遞歸
            for (t = T -1; t >= 1; t--)
            {
                for (i = 1; i <= phmm->N; i++)
                {
                    sum = 0.0;
                    for (j = 1; j <= phmm->N; j++)
                          sum += phmm->A[i][j] ＊ (phmm->B[j][O[t+1]])＊beta[t+1][j];
                    beta[t][i] = sum;
                }
            }
            //終止
            ＊pprob = 0.0;
            for (i = 1; i <= phmm->N; i++)
                ＊pprob += beta[1][i];
        }
        void BackwardWithScale(HMM ＊phmm, int T, int ＊O, double ＊＊beta, double ＊scale,double ＊pprob)
        {
            int i, j;
            int t;
            double sum;
            //初始化
            for (i = 1; i <= phmm->N; i++)
                beta[T][i] = 1.0/scale[T];
            //遞歸
            for (t = T -1; t >= 1; t--)
            {
                for (i = 1; i <= phmm->N; i++)
                {
                    sum = 0.0;
                    for (j = 1; j <= phmm->N; j++)
                        sum += phmm->A[i][j] ＊ (phmm->B[j][O[t+1]])＊beta[t+1][j];
                    beta[t][i] = sum/scale[t];
                }
            }
        }

Viterbi算法的代碼如代碼2-5所示。

代碼2-5 Viterbi算法代碼

        #define VITHUGE  100000000000.0
        void Viterbi(HMM ＊phmm, int T, int ＊O, double ＊＊delta, int ＊＊psi, int ＊q, double ＊pprob)
        {
            int i, j;
            int t;
            int maxvalind;
            double  maxval, val;
            //初始化
            for (i = 1; i <= phmm->N; i++)
            {
                delta[1][i] = phmm->pi[i] ＊ (phmm->B[i][O[1]]);
                psi[1][i] = 0;
            }
            //遞歸
            for (t = 2; t <= T; t++)
            {
                for (j = 1; j <= phmm->N; j++)
                {
                    maxval = 0.0;
                    maxvalind = 1;
                    for (i = 1; i <= phmm->N; i++)
                    {
                          val = delta[t-1][i]＊(phmm->A[i][j]);
                          if (val > maxval) {
                              maxval = val;
                              maxvalind = i;
                        }
                    }
                    delta[t][j] = maxval＊(phmm->B[j][O[t]]);
                    psi[t][j] = maxvalind;
                }
            }
            //終止
            ＊pprob = 0.0;
            q[T] = 1;
            for (i = 1; i <= phmm->N; i++)
            {
                if (delta[T][i] > ＊pprob)
                {
                    ＊pprob = delta[T][i];
                    q[T] = i;
                }
            }
            for (t = T -1; t >= 1; t--)
                q[t] = psi[t+1][q[t+1]];
      }
      void ViterbiLog(HMM ＊phmm, int T, int ＊O, double ＊＊delta, int ＊＊psi, int ＊q, double
            ＊pprob)
      {
            int i, j;
            int t;
            int maxvalind;
            double  maxval, val;
            double  ＊＊biot;
            //預處理
            for (i = 1; i <= phmm->N; i++)
                phmm->pi[i] = log(phmm->pi[i]);
            for (i = 1; i <= phmm->N; i++)
                for (j = 1; j <= phmm->N; j++)
                {
                    phmm->A[i][j] = log(phmm->A[i][j]);
                }
            biot = dmatrix(1, phmm->N, 1, T);
            for (i = 1; i <= phmm->N; i++)
                for (t = 1; t <= T; t++)
                {
                    biot[i][t] = log(phmm->B[i][O[t]]);
                }
            //初始化
            for (i = 1; i <= phmm->N; i++)
            {
                delta[1][i] = phmm->pi[i] + biot[i][1];
                psi[1][i] = 0;
            }
            //遞歸
            for (t = 2; t <= T; t++)
            {
                for (j = 1; j <= phmm->N; j++)
                {
                    maxval = -VITHUGE;
                    maxvalind = 1;
                    for (i = 1; i <= phmm->N; i++)
                    {
                        val = delta[t-1][i] + (phmm->A[i][j]);
                        if (val > maxval)
                        {
                              maxval = val;
                              maxvalind = i;
                        }
                    }
                    delta[t][j] = maxval + biot[j][t];
                    psi[t][j] = maxvalind;
                }
            }
            //終止
            ＊pprob = -VITHUGE;
            q[T] = 1;
            for (i = 1; i <= phmm->N; i++)
            {
                if (delta[T][i] > ＊pprob)
                {
                    ＊pprob = delta[T][i];
                    q[T] = i;
                }
            }
            //回朔
            for (t = T -1; t >= 1; t--)
                q[t] = psi[t+1][q[t+1]];
        }

Baum-Welch算法的代碼如代碼2-6所示。

代碼2-6 Baum-Welch算法代碼

        #define DELTA 0.001
        void BaumWelch(HMM ＊phmm, int T, int ＊O, double ＊＊alpha, double ＊＊beta, double＊＊gamma, int ＊pniter, double ＊plogprobinit, double ＊plogprobfinal)
        {
                int i, j, k;
                int t, l = 0;
                double  logprobf, logprobb;
                double  numeratorA, denominatorA;
                double  numeratorB, denominatorB;
                double ＊＊＊xi, ＊scale;
                double delta, deltaprev, logprobprev;
                deltaprev = 10e-70;
                xi = AllocXi(T, phmm->N);
                scale = dvector(1, T);
                ForwardWithScale(phmm, T, O, alpha, scale, &logprobf);
                ＊plogprobinit = logprobf;
                BackwardWithScale(phmm, T, O, beta, scale, &logprobb);
                ComputeGamma(phmm, T, alpha, beta, gamma);
                ComputeXi(phmm, T, O, alpha, beta, xi);
                logprobprev = logprobf;
                do{
         //重新估計t=1時，狀態為i的頻率
         for (i = 1; i <= phmm->N; i++)
phmm->pi[i] = .001 + .999＊gamma[1][i];
         //重新估計轉移矩陣和觀察矩陣
         for (i = 1; i <= phmm->N; i++)
         {
denominatorA = 0.0;
for (t = 1; t <= T -1; t++)
denominatorA += gamma[t][i];
for (j = 1; j <= phmm->N; j++)
{
numeratorA = 0.0;
for (t = 1; t <= T -1; t++)
numeratorA += xi[t][i][j];
phmm->A[i][j] = .001 + .999＊numeratorA/denominatorA;
}
denominatorB = denominatorA + gamma[T][i];
for (k = 1; k <= phmm->M; k++)
{
numeratorB = 0.0;
for (t = 1; t <= T; t++)
{
if (O[t] == k)
numeratorB += gamma[t][i];
                        }
                        phmm->B[i][k] = .001 +.999＊numeratorB/denominatorB;
                    }
                }
                ForwardWithScale(phmm, T, O, alpha, scale, &logprobf);
                BackwardWithScale(phmm, T, O, beta, scale, &logprobb);
                ComputeGamma(phmm, T, alpha, beta, gamma);
                ComputeXi(phmm, T, O, alpha, beta, xi);
                //計算兩次直接的概率差
                delta = logprobf - logprobprev;
                logprobprev = logprobf;
                l++;
            }
            while (delta > DELTA);
            ＊pniter = l;
            ＊plogprobfinal = logprobf; /＊ log P(O|estimated model) ＊/
            FreeXi(xi, T, phmm->N);
            free_dvector(scale, 1, T);
      }
      void ComputeGamma(HMM ＊phmm, int T, double ＊＊alpha, double ＊＊beta, double ＊＊gamma)
      {
            int i, j;
            int t;
            double  denominator;
            for (t = 1; t <= T; t++)
            {
                denominator = 0.0;
                for (j = 1; j <= phmm->N; j++)
                {
                    gamma[t][j] = alpha[t][j]＊beta[t][j];
                    denominator += gamma[t][j];
                }
                for (i = 1; i <= phmm->N; i++)
                    gamma[t][i] = gamma[t][i]/denominator;
            }
      }
      void ComputeXi(HMM＊ phmm, int T, int ＊O, double ＊＊alpha, double ＊＊beta, double ＊＊＊xi)
      {
            int i, j;
            int t;
            double sum;
            for (t = 1; t <= T -1; t++) {
                sum = 0.0;
                for (i = 1; i <= phmm->N; i++)
                    for (j = 1; j <= phmm->N; j++)
                    {
                        xi[t][i][j] = alpha[t][i]＊beta[t+1][j]
                              ＊(phmm->A[i][j])
                              ＊(phmm->B[j][O[t+1]]);
                        sum += xi[t][i][j];
                    }
                for (i = 1; i <= phmm->N; i++)
                    for (j = 1; j <= phmm->N; j++)
                        xi[t][i][j]  /= sum;
            }
        }
        double ＊＊＊ AllocXi(int T, int N)
        {
            int t;
            double ＊＊＊xi;
            xi = (double ＊＊＊) malloc(T＊sizeof(double ＊＊));
            xi --;
            for (t = 1; t <= T; t++)
                xi[t] = dmatrix(1, N, 1, N);
            return xi;
        }
        void FreeXi(double ＊＊＊ xi, int T, int N)
        {
            int t;
            for (t = 1; t <= T; t++)
                free_dmatrix(xi[t], 1, N, 1, N);
            xi ++;
            free(xi);
        }

官术网_书友最值得收藏!

Visual C++數字圖像模式識別典型案例詳解

2.2.2 隱馬爾可夫模型