1時間ごとに更新!Amazon.co.jpで今売れている本トップ100   掲示板ランキング



掲示板利用宣言

 次のフォームをすべてチェックしてからご利用ください。

 私は

 題名と投稿者名は具体的に書きます。
 課題の丸投げはしません。
 ソースの添付は「HTML変換ツール」で字下げします。
 返信の引用は最小限にします。
 環境(OSとコンパイラ)や症状は具体的に詳しく書きます。
 返信の付いた投稿は削除しません。
 マルチポスト(多重投稿)はしません。

掲示板1

管理者用メニュー    ツリーに戻る    携帯用URL    ホームページ    記事検索    ログ    タグ一覧

No.5692

Viterbiアルゴリズムの実装
投稿者---ジョン(2006/05/06 21:52:05)


隠れマルコフモデル(以下、HMM)を構築し、
Viterbiアルゴリズムを用いて
そのHMMにおける任意の文字列(アミノ酸配列)の最大確率を求めるために
以下のソースを作成したのですが、何度実行しても、
実際の最大確率よりも低い値が出てしまいます。

その原因がどこにあるのか、検討がつきません。
誤っている、またはこうした方がより簡潔で分かりやすいという点があれば、是非、ご教授下さい。

OSは、Windows Xp Home edition, コンパイラはgccを使っています。

#include <stdio.h>
#include <stdlib.h>  /* 擬似乱数操作関数を扱うヘッダファイルのインクルード */
#include <sys/time.h> /* PCより時間を取得する関数を扱うヘッダファイルをインクルード */
#include <math.h>  /* 数学関数を扱うヘッダファイルのインクルード */

/**************/
/* 定数の定義 */
/**************/
#define K 13      /* 状態の数 */
#define Z 5    /* 配列の長さ */
#define NUM_A 20    /* アミノ酸の種類の数 */
#define TYPE_S 0    /* 開始状態 */
#define TYPE_N 1    /* 通常のアミノ酸を出力する状態 */
#define LOG0 -1000.0  /* 無限小の意味でLOG0を定義 */

/**************************/
/* 関数のプロトタイプ宣言 */
/**************************/
double viterbi(int l, int i);
int a_to_n(char x);
double rand01();

/************************/
/* グローバル変数の宣言 */
/************************/

/* Viterbi アルゴリズムを適用するアミノ酸配列 */
static char x[Z];

/* 状態のタイプ */
static int type[K]={0};

/* 遷移確率の行列 */
static double t[K][K]={0.0};

/* 記号出力確率の行列 */
static double e[K][NUM_A]={0.0};

/************/
/* maim関数 */
/************/
int main(void){
  int l, i;
  int a, b;
  double sum = 0.0;
  struct timeval tv;
  gettimeofday(&tv,NULL);
  srand(tv.tv_usec);
  
  /* 状態のタイプ */
  for(a = 1; a < K; a++){
    type[a] = TYPE_N;
  }
  
  /* 配列の入力 */
  printf("Please, input sequence(Length = %d).\n sequence: ", Z);
  scanf("%s", x);
  
  /* 遷移確率の初期条件 */ 
  for(a = 0; a < K-1; a = a+2){
      t[a][a+1] = rand01();
      t[a][a+2] = 1 - t[a][a+1];
      t[a+1][a+1] = rand01();
      t[a+1][a+2] = 1 - t[a+1][a+1];
  }
  for(a = 0; a < K; a++){
    t[K-1][a] = 0;
  }
  
  /* 出力確率の初期条件 */
  for(a = 1; a < K-1; a++){
    for(b = 0; b < NUM_A; b++){
      e[a][b] = rand01();
      sum = sum + e[a][b];
    }
    for(b = 0; b < NUM_A; b++){
      e[a][b] /= sum;
    }
    sum = 0.0;
  }
  printf("\n");
  
  /* 遷移確率の出力*/
  printf("transition probability:\n");
  printf("|");
  for(a = 0; a < K; a++){
    printf("   %d  |",a);
  }
  printf("\n");
  
  for(a = 0; a < K; a++){
    printf("-------");
  }
  printf("\n");
  
  for(a = 0; a < K; a++){
    for(b = 0; b < K; b++){
      printf("  %.2f ", t[a][b]);
    }
    printf("| %d\n", a);
  }
  printf("\n");
  
  /* 出力確率の出力 */
  printf("output probability:\n");
  printf("|   A  |   C  |   D  |   E  |   F  |   G  |   H  |   I  |   K  |   L  |   M  |   N  |   P  |   Q  |   R  |   S  |   T  |   V  |   W  |   Y  |\n");
  printf("------------------------------------------------------------------------------------------------------------------------------------------------\n");
  for(a = 0; a < K; a++){
    for(b = 0; b < NUM_A ; b++){
      printf("  %.2f ", e[a][b]);
    }
    printf("| %d\n", a);
  }
  printf("\n");
  
  /* 状態の入力 */
  printf("Please, input state(0〜%d).\n state = ", K-1);
  scanf("%d",&l);
  printf("\n");
  
  /* 配列の何番目かを入力 */
  printf("Please, input ordinal number(1〜%d).\n ordinal number = ",Z);
  scanf("%d",&i);
  
  /* 最大確率の出力 */
  printf("Viterbi(l,i) = %.2f\n",viterbi(l,i));
  printf("\n");
}

/************************/
/* Viterbi アルゴリズム */
/************************/
double viterbi(int l, int i){
  int m, m_max;  /* 状態lの一つ前の状態 */
  double p, p_max;  /* 状態kに至るまでの確率 */  
  
  switch(type[l]){  /* 状態lのタイプを判定 */
    case TYPE_S:  /* 開始状態の場合 */
      if(i <= 0){ /* 開始状態はアミノ酸が出力されないので確率1,つまりその対数は0 */
        return 0.0;
      }
      else{  /* 開始状態はアミノ酸は出力されないのに出力される確率は0 */
        return LOG0;
      }
      break;
      
    case TYPE_N:  /* 通常の状態の場合 */
      if(i <= 0){
        return LOG0;
      }
      else{
        for(p_max = LOG0 -1, m = 0; m <=l; m++){  /* 最大の確率を与える状態 k_max およびその確率 p_max を探索 */
          if(t[m][l] > 0.0){
            p = viterbi(m, i-1) + log(t[m][l]);
          }
          else{
            p = LOG0;
          }
          if(p > p_max){p_max = p; m_max = m;}  /* p_maxの更新 */
        }
      }
      return p_max + log(e[l][a_to_n(x[i-1])]);
      break;
  }
}

/************************/
/* アミノ酸を数字に変換 */
/************************/
int a_to_n(char x){
  switch(x){
    case'A': return 0;
    case'C': return 1;
    case'D': return 2;
    case'E': return 3;
    case'F': return 4;
    case'G': return 5;
    case'H': return 6;
    case'I': return 7;
    case'K': return 8;
    case'L': return 9;
    case'M': return 10;
    case'N': return 11;
    case'P': return 12;
    case'Q': return 13;
    case'R': return 14;
    case'S': return 15;
    case'T': return 16;
    case'V': return 17;
    case'W': return 18;
    case'Y': return 19;
  }
}

/************************/
/* 0〜1の擬似乱数の生成 */
/************************/
double rand01()
{
  return (double)rand()/RAND_MAX;
}




この投稿にコメントする

削除パスワード

発言に関する情報 題名 投稿番号 投稿者名 投稿日時
<子記事> Viterbiアルゴリズムの実装 5693 ジョン 2006/05/06 21:52:50
<子記事> Re:Viterbiアルゴリズムの実装 5694 ぽへぇ 2006/05/07 07:34:26


No.5693

Viterbiアルゴリズムの実装
投稿者---ジョン(2006/05/06 21:52:50)


以下は実行例です。
Please, input sequence(Length = 5).
 sequence: AAAAA

transition probability:
|   0  |   1  |   2  |   3  |   4  |   5  |   6  |   7  |   8  |   9  |   10  |   11  |   12  |
-------------------------------------------------------------------------------------------
  0.00   0.98   0.02   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00 | 0
  0.00   0.58   0.42   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00 | 1
  0.00   0.00   0.00   0.82   0.18   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00 | 2
  0.00   0.00   0.00   0.83   0.17   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00 | 3
  0.00   0.00   0.00   0.00   0.00   0.64   0.36   0.00   0.00   0.00   0.00   0.00   0.00 | 4
  0.00   0.00   0.00   0.00   0.00   0.80   0.20   0.00   0.00   0.00   0.00   0.00   0.00 | 5
  0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.71   0.29   0.00   0.00   0.00   0.00 | 6
  0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.11   0.89   0.00   0.00   0.00   0.00 | 7
  0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.90   0.10   0.00   0.00 | 8
  0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.36   0.64   0.00   0.00 | 9
  0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.98   0.02 | 10
  0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.22   0.78 | 11
  0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00 | 12

output probability:
|   A  |   C  |   D  |   E  |   F  |   G  |   H  |   I  |   K  |   L  |   M  |   N  |   P  |   Q  |   R  |   S  |   T  |   V  |   W  |   Y  |
------------------------------------------------------------------------------------------------------------------------------------------------
  0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00 | 0
  0.08   0.01   0.03   0.07   0.02   0.01   0.09   0.05   0.05   0.08   0.08   0.03   0.08   0.07   0.00   0.08   0.05   0.06   0.06   0.02 | 1
  0.05   0.10   0.03   0.02   0.08   0.09   0.01   0.07   0.01   0.02   0.06   0.09   0.02   0.07   0.01   0.05   0.06   0.02   0.09   0.06 | 2
  0.04   0.04   0.07   0.06   0.07   0.08   0.09   0.02   0.07   0.01   0.08   0.03   0.06   0.05   0.07   0.00   0.06   0.03   0.01   0.04 | 3
  0.08   0.10   0.07   0.11   0.10   0.03   0.00   0.03   0.04   0.00   0.07   0.01   0.04   0.12   0.07   0.01   0.01   0.00   0.03   0.08 | 4
  0.07   0.07   0.06   0.01   0.05   0.01   0.10   0.06   0.06   0.06   0.06   0.03   0.03   0.00   0.00   0.09   0.08   0.03   0.05   0.08 | 5
  0.00   0.04   0.02   0.04   0.00   0.01   0.10   0.03   0.03   0.02   0.11   0.04   0.03   0.03   0.10   0.07   0.00   0.11   0.12   0.11 | 6
  0.07   0.08   0.04   0.02   0.02   0.00   0.08   0.08   0.07   0.04   0.08   0.06   0.08   0.07   0.04   0.08   0.05   0.00   0.04   0.00 | 7
  0.08   0.07   0.04   0.05   0.05   0.08   0.03   0.07   0.04   0.02   0.01   0.04   0.04   0.01   0.06   0.05   0.04   0.07   0.07   0.08 | 8
  0.00   0.08   0.04   0.12   0.06   0.07   0.06   0.04   0.05   0.09   0.02   0.03   0.00   0.07   0.03   0.06   0.06   0.02   0.03   0.07 | 9
  0.05   0.03   0.09   0.03   0.04   0.04   0.05   0.06   0.08   0.01   0.07   0.00   0.06   0.08   0.06   0.07   0.06   0.00   0.07   0.06 | 10
  0.04   0.07   0.06   0.08   0.04   0.06   0.08   0.11   0.02   0.02   0.00   0.10   0.03   0.08   0.02   0.05   0.09   0.00   0.03   0.02 | 11
  0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00 | 12

Please, input state(0〜12).
 state = 1

Please, input ordinal number(1〜5).
 ordinal number = 1
Viterbi(l,i) = -2.58



この投稿にコメントする

削除パスワード

No.5694

Re:Viterbiアルゴリズムの実装
投稿者---ぽへぇ(2006/05/07 07:34:26)


>実際の最大確率よりも低い値が出てしまいます。
期待している最大確率ってどのくらいなのですか?
また、それ(期待している最大確率)はどのような方法で
求めたのですか?

その方法とプログラムが異なる点はどこでしょうか?
異なる点を見つけ出すにはどうしたらよいのでしょうか?



この投稿にコメントする

削除パスワード

No.5695

Re:Viterbiアルゴリズムの実装
投稿者---ジョン(2006/05/07 12:03:57)


ぽへぇさん、返信ありがとうございます。

>期待している最大確率ってどのくらいなのですか?

No.5693 の場合、-1.11 のはずが -2.58 と表示されています。



>また、それ(期待している最大確率)はどのような方法で
>求めたのですか?

このプログラム内の HMM は状態 0 から出発し、
状態 12 まで遷移して、5 アミノ酸以上の配列の
最大確率を求められるようにしたつもりです。
その確率は状態kから状態lへの遷移確率 t[k][l] と
状態 l でのアミノ酸の出力確率 e[l] から求められます。

No.5693 の場合、状態 1 で 1 番目のアミノ酸
(入力された配列が AAAAA なので A)が出力される最大確率は
状態 0 から状態 1 への遷移で 0.98 (これは transition probability の 1 行 2 列目)、
状態 1 で A を出力する確率は 0.08 (これは output probability の 2 行 1 列目)です。
したがって、
viterbi(1, 1) = log(0.98 * 0.08) = -1.11
となります。



>その方法とプログラムが異なる点はどこでしょうか?

まさにその点をお聞きしたいのです。
計算過程を手書きで確認してからソースを書いたのですが、
異なる点がどこなのかがわからず、皆さんに質問させて頂きました。



>異なる点を見つけ出すにはどうしたらよいのでしょうか?

関数 viterbi(int l, int i) 内で戻り値である
viterbi(l, i)を出力させ、
v(0, 0), v(a, 0) (a > 0), v(0, b) (b > 0)
は正常な値を示していたのですが、これら以外は異なる値が
出力されてしまいます。


この投稿にコメントする

削除パスワード

No.5696

Re:Viterbiアルゴリズムの実装
投稿者---ぽへぇ(2006/05/07 12:25:10)



関数 log() と log10() の違いを調べると幸せになれるかも。



この投稿にコメントする

削除パスワード

No.5697

Re:Viterbiアルゴリズムの実装
投稿者---ジョン(2006/05/07 12:37:16)


自然対数で計算していたなんて…
ぽへぇさん、本当にありがとうございました。


この投稿にコメントする

削除パスワード

管理者用メニュー    ツリーに戻る    携帯用URL    ホームページ    記事検索    ログ    タグ一覧





掲示板提供:(有)リアル・インテグリティ