大資料之道 HMM系列

標籤：維特比演算法   hmm   viterbi   自然語言處理   

一：HMM解碼問題

（1）給定一個觀察序列O=O1O2...OT,和模型μ=（A,B,π），如何快速有效地選擇在一定意義下“最優”的狀態序列Q=q1q2...qT，使該狀態最好地解釋觀察序列。

（2）最可能的隱藏狀態序列(Finding most probable sequence of hidden states)；對於一個特殊的隱馬爾科夫模型(HMM)及一個相應的觀察序列，我們常常希望能找到產生此序列最可能的隱藏狀態序列。

二：執行個體篇

（1）假設連續觀察3天的海藻濕度為(Dry,Damp,Soggy),求這三天最可能的天氣情況。天氣只有三類(Sunny,Cloudy,Rainy)，而且海藻濕度和天氣有一定的關係。

已知：
1. 隱藏的狀態：Sunny,Cloudy, Rainy;海藻濕度有四類{Dry,Dryish,Damp,Soggy }
2. 觀察狀態序列：{Dry, Damp, Soggy };
3. 初始狀態序列：Sunny(0.63),Cloudy(0.17),Rainy(0.20);

4. 狀態轉移矩陣

   Sunny        Cloudy   Rainy
Sunny 0.5         0.375  0.125
Cloudy 0.25 0.125 0.625
Rainy 0.25 0.375 0.375

Cloudy(昨天)->Sunny(今天)的機率是0.25;

Sunny(昨天)->Rainy(今天)的機率是0.125.

5. 混淆矩陣（海藻濕度與天氣的相關性）：

         Dry          Dryish DampSoggy
Sunny 0.6         0.2         0.150.05
Cloudy 0.25 0.25 0.25 0.25
Rainy 0.05 0.10 0.35 0.50

（2）計算方法：

由HMM可知，Day2的天氣僅取決於Day1；Day3的天氣又只取決於Day2的天氣。

step1：Day1由於是初始狀態，我們分別求

P(Day1-Sunny)=0.63*0.6;
P(Day1-Cloudy)=0.17*0.25;
P(Day1-Rain)=0.20*0.05;
Choose max{ P(Day1-Sunny) , P(Day1-Cloudy),P(Day1-Rainy)}, 得到P(Day1-Sunny)最大，得出第1天Sunny的機率最大。

step2: Day2的天氣又取決於Day1的天氣狀況，同時也受Day2觀察的海藻情況影響。

P(Day2-Sunny)= max{ P(Day1-Sunny)*0.5, P(Day1-Cloudy)*0.25,  P(Day1-Rainy)*0.25} *0.15;
P(Day2-Cloudy)= max{ P(Day1-Sunny)*0.375,  P(Day1-Cloudy)*0.125, P(Day1-Rainy)*0.625} *0.25;
P(Day2-Rainy)= max{ P(Day1-Sunny)*0.125,  P(Day1-Cloudy)*0.625 , P(Day1-Rainy)*0.375} *0.35;
Choosemax{ P(Day2-Sunny) , P(Day2-Cloudy), P(Day2-Rainy)},得到P(Day2-Rainy)最大，得出第2天Rainy的機率最大。
故{Sunny,Rainy}是前兩天最大可能的天氣序列。

step3: Day3的天氣又取決於Day2的天氣狀況，同時也受Day3觀察的海藻情況影響。

P(Day3-Sunny)= max{ P(Day2-Sunny)*0.5, P(Day2-Cloudy)*0.25,  P(Day2-Rainy)*0.25} *0.05;
P(Day3-Cloudy)= max{ P(Day2-Sunny)*0.375,  P(Day2-Cloudy)*0.125, P(Day2-Rainy)*0.625} *0.25;
P(Day3-Rainy)= max{ P(Day2-Sunny)*0.125,  P(Day2-Cloudy)*0.625, P(Day2-Rainy)*0.375} *0. 05;
Choosemax{ P(Day3-Sunny) , P(Day3-Cloudy), P(Day3-Rainy)},得到 P(Day3-Rainy)最大，得出第3天Rainy的機率最大。

(3)代碼實現

#include <stdio.h>#include <stdlib.h>#include <string.h>/* 題目描述：給定一個觀察序列O=O1O2...OT,和模型μ=（A,B,π），如何快速有效地選擇在一定意義下“最優”的狀態序列Q=q1q2...qT，使該狀態最好地解釋觀察序列。下面a是狀態轉移矩陣； b是發射矩陣，list 是觀察序列；*/const int HID_STATUSES = 3;  // 隱藏的狀態個數const int VIS_ST = 4;        // 觀察到的實際所有可能的情況（狀態）個數const int N = 3;             // 觀察的天數（序列長度）int main(){    float a[HID_STATUSES][HID_STATUSES] = {{0.5,0.375,0.125},{0.25,0.125,0.625},{0.25,0.375,0.375}};//狀態轉移矩陣    float b[HID_STATUSES][VIS_ST] = {{0.6,0.2,0.15,0.05},{0.25,0.25,0.25,0.25,},{0.05,0.10,0.35,0.50,}};//b是發射矩陣    float ini[HID_STATUSES] = {0.63,0.17,0.20};// 初始的    float result[N][HID_STATUSES];    int list[N] = {0,2,3};// 測試用的觀察到的狀態序列,4天的觀測值；由於就有2中觀察狀態，所以用    // 0 和 1 區分，實際中應該是input的，以及觀察的天數也是input的或者read file    int max[N][HID_STATUSES];    float tmp;    int i,j,k;    //step1:Initialization，其中{0.2,0.4,0.4} 是π    for(i=0;i<HID_STATUSES;++i)    {        result[0][i] = ini[i]*b[i][list[0]];    }    int count = 1, max_node;// count 是list的下標，從1開始，第一天是step1已經初始化    float max_v;    //step2:歸納運算 ，i 從1開始，應為step1已經初始化了    for (i=1; i<N; i++)    {        for(j=0; j<HID_STATUSES; j++)        {            tmp = result[i-1][0] * a[0][j] * b[j][list[count]];            max[i][j] = 0;            for(k=1; k<HID_STATUSES; k++)            {                if(result[i-1][k] * a[k][j] * b[j][list[count]] > tmp)                {                    tmp = result[i-1][k] * a[k][j]* b[j][list[count]];                    max[i][j] = k;                }               result[i][j] = tmp;            }            max_v = result[N-1][0];            max_node = 0;            for (k=1; k<HID_STATUSES; k++)            {                if(result[N-1][k] > max_v)                {                    max_v = result[N-1][k];                    max_node = k;                }            }        }        count += 1;    }    //step3:終結   for (i=0; i<N; i++)    {        for(j=0; j<HID_STATUSES; j++)        {            printf("%d %d     %f\n",i+1,j+1,result[i][j]);        }    }    printf("Pmax= %f\n", max_v);    printf("step%d: %d   \n",N,max_node+1);    //step4:回溯    for(k=N-1; k>=1; k--)    {        printf("step%d: %d  \n",k, max[k][max_node]+1);        max_node = max[k][max_node];    }    return 0;}

（4）運行結果：

大資料之道 HMM系列