Crfと素性テンプレート
View more presentations from uchumik
TokyoNLP で話す内容を貼っておきます。
きっと後で訂正が入るんだろうな。。
2011年6月25日土曜日
2011年6月12日日曜日
実数素性テンプレートの使える CRF
この間記事で書いた実数素性テンプレートを使える CRF を作ってみました。
google code
面倒くさかったのでだいぶコピペの目立つコードになってますが、その辺は勘弁。
以下は動作確認のためだけのコード。
下が動かしたときの出力。
実数素性を使えそうな公開データとかあればそれで試してみたいですけど、とりあえず今日の所はここまで。
google code
面倒くさかったのでだいぶコピペの目立つコードになってますが、その辺は勘弁。
- テンプレートが文字列を扱うバイナリな値を取る場合（%x マクロを使うもの）ならば、素性関数 f_k(y_{i-1},y_i,X) は 1 or 0 を取る。
- テンプレートが実数を扱う場合（%r や SUM, MAX, Value, etc. ）ならば、素性関数 f_k(y_{i-1},y_i,X) はテンプレートに書かれた演算結果の値 or 0 を取る。
以下は動作確認のためだけのコード。
# include "rtcrflearn.hpp"
# include "rtcrftagger.hpp"
using namespace RtCrf;
using namespace std;
int main(int argc, char **argv)
{
Crflearn learner(*(argv+1),*(argv+2),1000000);
learner.init();
learner.learn(5,0);
learner.save("test.model");
Crftagger tagger(*(argv+1),1000000);
tagger.read("test.model");
tagger.tagging(*(argv+2));
return 0;
}
下が動かしたときの出力。
# ./test template ../data/train.txt
labels: 2
bound: 3
ufeatures: 17
bfeatures: 1
instance: 1
uparameters: 34
bparameters: 8
model parameters: 42
epoch: 0 err:0.500000(2/4)
epoch: 1 err:0.000000(0/4)
epoch: 2 err:0.000000(0/4)
epoch: 3 err:0.000000(0/4)
epoch: 4 err:0.000000(0/4)
a b c 1.0 -1.0 L_A L_A
d e f -1.0 1.0 L_B L_B
p e f -0.5 0.1 L_B L_B
g b c 0.3 -0.1 L_A L_A
実数素性を使えそうな公開データとかあればそれで試してみたいですけど、とりあえず今日の所はここまで。
Tagger クラスの出力に予測時のスコアを出すようにすれば、１段目の CRF の出力を２段目の CRF の学習素性に使うというようなスタッキングも出来るようになるかなと思う次第。
そのうち CNF とかと一緒にまとめて API を変えとこうかな。
2011年6月8日水曜日
実数素性テンプレートを作ろう。
まだ生きてます。
ここしばらくはちと研究のまねごとなどを会社でしている事もあって、あんまりブログに書ける事も無い日々です。
とりあえず対外的に問題ない範囲だと、去年の11月にACMLという会議でポスター発表をしてみました。
ポスターはアブストの提出だけだったのですが、参考文献が後で見直したら壊れていたという大失敗。
後、同じ11月にNL研でクエリの訂正手法についての発表をしてきました。
広島のお好み焼きはおいしかったです。
あと、 NetWalker を活用するべくいくつか実装したものもあるのですが、そちらはまだしばらく非公開。
とりあえず生存報告はこんな所。
最近忙しいですが何となく以前よりは楽しく過ごせている気がします。
で、久しぶりに書いた記事がこれで終わりというのもなんなので、近々 TokyoNLP で発表しようかなと思っている CRF の素性テンプレート周りの話に関連して少々追記。
素性テンプレートが使える CRF というと、CRF++ や crfsgd がありますが、基本的にこれらの素性テンプレートは文字列を素性とするタイプで、素性関数は f_k(y,x) = 1 or 0 というようにバイナリの値を取るようになっています。
これはある素性関数は特定の y と x のペアが観測された時に 1 を返して、それ以外のときは 0 を取るという意味で、素性がテキストの時は特に問題になる事はありません。
が、最近はクエリログとかを扱う機会も多くなってきた都合上、出来れば実数値も扱いたいなーと思う事がしばしばです。
例えば、クエリに含まれるトークンにタグ付けをする際には、トークンの尤度を言語モデルで求めて、それを素性として扱いたいとかがあると思うのですが、それは現状だと数値の範囲を表す様なテキストに変更して素性として扱うしかなさそうです。
ですが、この方法だとスコアの範囲ごとに素性が分かれるので、やっぱり元々の言語モデルで求めた尤度を直接素性として扱うのとは違いそうです。
なので、実数の扱える素性テンプレートを作ってみました。
http://code.google.com/p/cnf/source/browse/#hg%2Ftoys%2Frtcrf%2Fsrc
動作としては、 f_k(y,テキスト素性) = 1 or 0、 f_k(y,実数素性) = 素性の実数値 or 0 となります。
もちろん、素性テンプレートというからには観測された素性の組み合わせを扱えないと格好が悪いので、実数素性でも組み合わせ（素性どうしの各種演算等）が扱えるようにしました。
で、ついでに素性テンプレートが今まで学習器とセットで独立性が低かった点も見直してみました。
使える演算としては、
をひとまず用意。
上の例だと x は変数に見えますが、別の演算を x の所に入れる事も出来るので、
みたいな事も出来ます。
とりあえず実際に適当に書いてみたテンプレートの例が以下です。
で、入力に使ったこれまた適当データが次。
これに上のテンプレートを使って素性を展開すると、次のようになります。
実装 は c++ で、 bison、 flex を使って作っています。
tsubosaka さんからはこのくらいの仕様なら bison、 flex なんて使わないでいいんじゃないの？というありがたいお言葉を頂いたのですが、残念ながら僕は 1 からパーザを書いた経験が無く、正直自信が無いためこれらを使っています。
ちなみに flex のバージョンが 2.5.35 、bison が 2.3 の環境で作っています。
bison は 2.1 以降なら c++ に対応しているので問題ないと思いますが、 flex の方は 2.5.35 だと yyFlexLexer::yywrap の実装がデフォルトで無いため、yywrap の実装を適当にソースコードの中に書いています。
ここしばらくはちと研究のまねごとなどを会社でしている事もあって、あんまりブログに書ける事も無い日々です。
とりあえず対外的に問題ない範囲だと、去年の11月にACMLという会議でポスター発表をしてみました。
ポスターはアブストの提出だけだったのですが、参考文献が後で見直したら壊れていたという大失敗。
後、同じ11月にNL研でクエリの訂正手法についての発表をしてきました。
広島のお好み焼きはおいしかったです。
あと、 NetWalker を活用するべくいくつか実装したものもあるのですが、そちらはまだしばらく非公開。
とりあえず生存報告はこんな所。
最近忙しいですが何となく以前よりは楽しく過ごせている気がします。
で、久しぶりに書いた記事がこれで終わりというのもなんなので、近々 TokyoNLP で発表しようかなと思っている CRF の素性テンプレート周りの話に関連して少々追記。
素性テンプレートが使える CRF というと、CRF++ や crfsgd がありますが、基本的にこれらの素性テンプレートは文字列を素性とするタイプで、素性関数は f_k(y,x) = 1 or 0 というようにバイナリの値を取るようになっています。
これはある素性関数は特定の y と x のペアが観測された時に 1 を返して、それ以外のときは 0 を取るという意味で、素性がテキストの時は特に問題になる事はありません。
が、最近はクエリログとかを扱う機会も多くなってきた都合上、出来れば実数値も扱いたいなーと思う事がしばしばです。
例えば、クエリに含まれるトークンにタグ付けをする際には、トークンの尤度を言語モデルで求めて、それを素性として扱いたいとかがあると思うのですが、それは現状だと数値の範囲を表す様なテキストに変更して素性として扱うしかなさそうです。
ですが、この方法だとスコアの範囲ごとに素性が分かれるので、やっぱり元々の言語モデルで求めた尤度を直接素性として扱うのとは違いそうです。
なので、実数の扱える素性テンプレートを作ってみました。
http://code.google.com/p/cnf/source/browse/#hg%2Ftoys%2Frtcrf%2Fsrc
動作としては、 f_k(y,テキスト素性) = 1 or 0、 f_k(y,実数素性) = 素性の実数値 or 0 となります。
もちろん、素性テンプレートというからには観測された素性の組み合わせを扱えないと格好が悪いので、実数素性でも組み合わせ（素性どうしの各種演算等）が扱えるようにしました。
で、ついでに素性テンプレートが今まで学習器とセットで独立性が低かった点も見直してみました。
使える演算としては、
- ( x )
- x + x
- x - x
- x * x
- x / x
- x % x
- -x
- log x
- sqrt x
- exp x
- max ( x,... )
- min ( x,... )
- sum ( x,... )
- prod ( x,... )
- x > x
- x >= x
- x < x
- x <= x
をひとまず用意。
上の例だと x は変数に見えますが、別の演算を x の所に入れる事も出来るので、
(x_1+x_2)*x_3*max(max(...),min(...),x_4,log(x_5),...)
みたいな事も出来ます。
とりあえず実際に適当に書いてみたテンプレートの例が以下です。
# Unigram
U00:%x[-2,0]
U01:%x[-1,0]_%x[0,0]
U02:%x[-1,0]_%x[0,0]_%x[1,0]
U03:1.0-%r[0,3]
U04:%r[0,3]+%r[0,4]
U05:%r[0,3]-%r[0,4]
U06:%r[0,3]*-%r[0,4]
U07:%r[0,3]/-%r[0,4]
U08:%r[0,3]%-%r[0,4]
U09:SUM(%r[.,3])
U10:PROD(%r[.,4])
U11:MIN(%r[0,3],%r[0,4])
U12:MAX(%r[0,3],%r[0,4])
U13:MIN(%r[.,3])
U14:MAX(%r[.,4])
U15:MIN(%r[0,.])
U16:MAX(%r[0,.])
U17:MIN(MAX(%r[0,.],%r[.,3]),%r[.,4])
U18:%r[0,3]/MAX(%r[.,3])
U19:%r[0,3]/SUM(%r[.,3])
U20:MIN(%r[0,3]+%r[0,4],%r[0,3])
# Bigram
B
で、入力に使ったこれまた適当データが次。
a b c 1.0 -1.0 L_A
d e f -1.0 1.0 L_B
p e f -0.5 0.1 L_B
g b c 0.3 -0.1 L_A
これに上のテンプレートを使って素性を展開すると、次のようになります。
U00:_B-2:1 U01:_B-1_a:1 U02:_B-1_a_d:1 U05:2 U06:1 U07:1 U09:-0.2 U10:0.01 U11:-1 U12:1 U13:-1 U14:1 U15:-1 U16:1 U17:-1 U18:1 U19:-5 B:1
U00:_B-1:1 U01:a_d:1 U02:a_d_p:1 U03:2 U05:-2 U06:1 U07:1 U09:-0.2 U10:0.01 U11:-1 U12:1 U13:-1 U14:1 U15:-1 U16:1 U17:-1 U18:-1 U19:5 U20:-1 B:1
U00:a:1 U01:d_p:1 U02:d_p_g:1 U03:1.5 U04:-0.4 U05:-0.6 U06:0.05 U07:5 U08:-0.1 U09:-0.2 U10:0.01 U11:-0.5 U12:0.1 U13:-1 U14:1 U15:-0.5 U16:0.1 U17:-1 U18:-0.5 U19:2.5 U20:-0.5 B:1
U00:d:1 U01:p_g:1 U02:p_g__E+1:1 U03:0.7 U04:0.2 U05:0.4 U06:0.03 U07:3 U08:0.1 U09:-0.2 U10:0.01 U11:-0.1 U12:0.3 U13:-1 U14:1 U15:-0.1 U16:0.3 U17:-1 U18:0.3 U19:-1.5 U20:0.2 B:1
実装 は c++ で、 bison、 flex を使って作っています。
tsubosaka さんからはこのくらいの仕様なら bison、 flex なんて使わないでいいんじゃないの？というありがたいお言葉を頂いたのですが、残念ながら僕は 1 からパーザを書いた経験が無く、正直自信が無いためこれらを使っています。
ちなみに flex のバージョンが 2.5.35 、bison が 2.3 の環境で作っています。
bison は 2.1 以降なら c++ に対応しているので問題ないと思いますが、 flex の方は 2.5.35 だと yyFlexLexer::yywrap の実装がデフォルトで無いため、yywrap の実装を適当にソースコードの中に書いています。
これが 2.5.4 以降だと標準で実装されているので、多分コンパイル時に多重定義で怒られるんじゃないかなぁという気がします。
今週はきっとこれ使った CRF 実装してます。
今週はきっとこれ使った CRF 実装してます。
2010年8月5日木曜日
bigram feature について。
twitter で @takeda25 さんが指摘されていたのですが, CRF で f(y_{i-1},y_{i},x_{i}) という観測素性とラベル bigram の組みまで考慮した素性関数があまり使われないのは, 素性数が増えるからというよりは計算に時間が掛かるためだろうか? という事を私なりに考えてみました.
CRF の1事例に対するパラメータ推定に掛かるオーダーはラベル数L, 系列長Tの時に, O(L^2T) です.
forward-backward でラティス中の位置 i, ラベル j のノードの alpha を計算する際には以下の logsumexp の計算を行います.
この計算を全ての i, j について行います.
この時足し込んでいる cost が 観測素性+遷移素性のコストで, 観測素性 x に対してラベル bigram を考慮する場合, 遷移素性のコスト計算をする時に一緒に足す事になるため, 単純に観測された bigram 素性の数だけ加算回数が増えます.
ですが, logsumexp の呼び出し回数自体は変わりません. (最初この部分も定数倍になると勘違いしていました)
cost の加算部分だけが定数倍される程度なので, 計算時間にはそんなに影響が出ないのではないかと考えたのですが, 実は CNF だと結構影響が出ました...
CNF では cost の計算時に, logstic 関数による非線形な重み付けが必要になるため, exp の計算が入ります. exp の計算が素性数だけ増えるため, 試してみたら結構遅くなりました..
試してみたのは CoNLL2000 の学習データで, 使用したテンプレートは以下の通りです. 見事に全部 bigram テンプレートです.
これで学習を行ってみたところ, 次のようになりました.
* iter = 1
ちなみに, テンプレートを元々の
にすると, iter = 1 の学習時間は以下の通り.
だいたい bigram テンプレート数倍だけ遅くなりました..
exp の計算が bigram テンプレート数倍必要になるので, そこが影響したようです.
CRF の場合は logistic 関数の計算無しで加算のみになるため, ここまでの影響は出ないと考えられるので, bigram 素性を使ったとしてもそこまで遅くなるという事は無いのではないかと思います.
ちなみにオール bigram テンプレートを使った時の学習に必要なメモリは約 450 MB でした.
bigram 素性を使うとタグ付け時にトークンの文脈依存性を上手く学習できたりしそうですし, 意味の無い素性ではないと思うのですが, 何であまり使われていないんでしょうね.
ちなみに, CRFSuite では高速化の為にいろいろな工夫がされているそうで, その工夫の中にラベル bigram と観測素性の組み合わせは考慮しないという前提に基づいた手法があるのか, bigram 素性を扱うと高速化手法のどれかが使えなくなってしまうっぽい?
ソースコードを落としてきてちら見してみたのですが, 人のコードを読むのが苦手な私ではすぐには問題箇所が分からず, 「おー! コード超キレイ!」とそんな感想しか言えませんでした...
時間のある時にじっくり読んで勉強させてもらいます!
CRF の1事例に対するパラメータ推定に掛かるオーダーはラベル数L, 系列長Tの時に, O(L^2T) です.
forward-backward でラティス中の位置 i, ラベル j のノードの alpha を計算する際には以下の logsumexp の計算を行います.
for (k = 0; k < L; ++k)
alpha_{i,j} = logsumexp(alpha_{i,j}, alpha_{i-1,k}+cost)
この計算を全ての i, j について行います.
この時足し込んでいる cost が 観測素性+遷移素性のコストで, 観測素性 x に対してラベル bigram を考慮する場合, 遷移素性のコスト計算をする時に一緒に足す事になるため, 単純に観測された bigram 素性の数だけ加算回数が増えます.
ですが, logsumexp の呼び出し回数自体は変わりません. (最初この部分も定数倍になると勘違いしていました)
cost の加算部分だけが定数倍される程度なので, 計算時間にはそんなに影響が出ないのではないかと考えたのですが, 実は CNF だと結構影響が出ました...
CNF では cost の計算時に, logstic 関数による非線形な重み付けが必要になるため, exp の計算が入ります. exp の計算が素性数だけ増えるため, 試してみたら結構遅くなりました..
試してみたのは CoNLL2000 の学習データで, 使用したテンプレートは以下の通りです. 見事に全部 bigram テンプレートです.
# Bigram
B00:%x[-2,0,0]
B01:%x[-1,0,0]
B02:%x[0,0,0]
B03:%x[1,0,0]
B04:%x[2,0,0]
B05:%x[-1,0,0]/%x[0,0,0]
B06:%x[0,0,0]/%x[1,0,0]
B10:%x[-2,1,0]
B11:%x[-1,1,0]
B12:%x[0,1,0]
B13:%x[1,1,0]
B14:%x[2,1,0]
B15:%x[-2,1,0]/%x[-1,1,0]
B16:%x[-1,1,0]/%x[0,1,0]
B17:%x[0,1,0]/%x[1,1,0]
B18:%x[1,1,0]/%x[2,1,0]
B20:%x[-2,1,0]/%x[-1,1,0]/%x[0,1,0]
B21:%x[-1,1,0]/%x[0,1,0]/%x[1,1,0]
B22:%x[0,1,0]/%x[1,1,0]/%x[2,1,0]
B
これで学習を行ってみたところ, 次のようになりました.
* iter = 1
time ./src/learner --learner=cnf -t src/template -c data/conll2000/train.txt -s bigramtest.save --l2 --pool=1000000 --block=24 --lambda=2 --iter=1
labels: 22
features: 338552
bound: 3
ufeatures: 0
bfeatures: 76329
instance: 8936
gate functions: 20
parameters of gate functions: 31
uparameters: 0
bparameters: 40301712
model parameters: 40301743
epoch: 0 err:0.069070(14624/211727)
./src/learner --learner=cnf -t src/template -c data/conll2000/train.txt -s 1323.90s user 41.84s system 96% cpu 23:40.72 total
ちなみに, テンプレートを元々の
U00:%x[-2,0,0]
U01:%x[-1,0,0]
U02:%x[0,0,0]
U03:%x[1,0,0]
U04:%x[2,0,0]
U05:%x[-1,0,0]/%x[0,0,0]
U06:%x[0,0,0]/%x[1,0,0]
U10:%x[-2,1,0]
U11:%x[-1,1,0]
U12:%x[0,1,0]
U13:%x[1,1,0]
U14:%x[2,1,0]
U15:%x[-2,1,0]/%x[-1,1,0]
U16:%x[-1,1,0]/%x[0,1,0]
U17:%x[0,1,0]/%x[1,1,0]
U18:%x[1,1,0]/%x[2,1,0]
U20:%x[-2,1,0]/%x[-1,1,0]/%x[0,1,0]
U21:%x[-1,1,0]/%x[0,1,0]/%x[1,1,0]
U22:%x[0,1,0]/%x[1,1,0]/%x[2,1,0]
# Bigram
B
にすると, iter = 1 の学習時間は以下の通り.
time ./src/learner --learner=cnf -t src/template -c data/conll2000/train.txt -s bigramtest.save --l2 --pool=1000000 --block=24 --lambda=2 --iter=1
labels: 22
features: 338552
bound: 3
ufeatures: 76328
bfeatures: 1
instance: 8936
gate functions: 20
parameters of gate functions: 31
uparameters: 1679216
bparameters: 528
model parameters: 1679775
epoch: 0 err:0.057314(12135/211727)
./src/learner --learner=cnf -t src/template -c data/conll2000/train.txt -s 67.23s user 5.61s system 92% cpu 1:18.67 total
だいたい bigram テンプレート数倍だけ遅くなりました..
exp の計算が bigram テンプレート数倍必要になるので, そこが影響したようです.
CRF の場合は logistic 関数の計算無しで加算のみになるため, ここまでの影響は出ないと考えられるので, bigram 素性を使ったとしてもそこまで遅くなるという事は無いのではないかと思います.
ちなみにオール bigram テンプレートを使った時の学習に必要なメモリは約 450 MB でした.
bigram 素性を使うとタグ付け時にトークンの文脈依存性を上手く学習できたりしそうですし, 意味の無い素性ではないと思うのですが, 何であまり使われていないんでしょうね.
ちなみに, CRFSuite では高速化の為にいろいろな工夫がされているそうで, その工夫の中にラベル bigram と観測素性の組み合わせは考慮しないという前提に基づいた手法があるのか, bigram 素性を扱うと高速化手法のどれかが使えなくなってしまうっぽい?
ソースコードを落としてきてちら見してみたのですが, 人のコードを読むのが苦手な私ではすぐには問題箇所が分からず, 「おー! コード超キレイ!」とそんな感想しか言えませんでした...
時間のある時にじっくり読んで勉強させてもらいます!
2010年7月9日金曜日
a hybrid markov/semi-markov cnf
名前だけ知りつつもどんな物なのか知らなかったので、5月に semi-markov crf の論文 を読んでみました。
semi-crf は入力系列に対する最適なセグメンテーションを学習する学習器です。
crf が P(Y|X, W) を求めていたのに対して、semi-crf では P(S|X, W) を推定します。
文章で書くと、2つの違いは
* crf では入力系列 X に対して生成されうるすべてのラベル系列 Y' に対して、正解系列とそれ以外の系列を弁別するよう学習する
* semi-crf では入力系列 X に対して生成されうるすべてのセグメンテーション結果 S' に対して、正解のセグメンテーションとそれ以外のセグメンテーションを弁別するよう学習する
ということになります。
semi-crf では入力系列が与えられた際に、t_j , ... , u_j までを1つのセグメント s_j として、1つのラベルを割り当てます。
この時、t_j 、 u_j は常に以下を満たす必要があります。
1 <= t_j <= u_j <= |S| , t_{j+1} = u_j + 1
これは、セグメントの終了位置と次のセグメントの開始位置は一致するよということです。
で、crf で使われていた素性関数 f の代わりに segment feature function g ( g^k(y_j,y_{j-1},X,t_j,u_j) )を使います。
といっても、あまり大きく変わっているという事は無くて、crf の素性関数がトークン x_i とラベル y_i のペアに対して 1 or 0 を返していたのが、セグメント s_i とラベル y_i のペアに変わった程度です。
大きな違いは以上なのですが、パラメータ推定をする際のラティスの形が違います。
crf ではラティスはきれいな格子でしたが、semi-crf では各位置 i について、i で始まるセグメントと i で終了するセグメントのみ接続を許すため、結構複雑な形をしています。
ただ、これは形態素解析の物と同じなので 自然言語処理-基礎と応用- の1章、形態素解析の 1.1.5 節あたりを読むと分かりやすいと思います。
パラメータ推定は crf と同じ方法が使えます。
とりあえずここまで理解したところで、cnf で作ってみよう！ と考えて作り始めたのですが、いざ作ってみて失敗に気づきました。
semi-crf ではセグメント単位で同じラベルをふるため、ラベル付けは IOB とかじゃ無く、トークン列に対して AAA、BBB のようにつければ良いやと思っていたのですが、いざ動かす所で、正解の与え方ではまりました。
というのも、AA, AA, AA というラベル列を正解コーパスでつけたつもりでも、AAA, AAA、AAAAAA のどれとも区別が出来なくなってました。。
それで、方向を修正し、結局ラベル付けは IOB タグに限定する事にして（cnf ではそんな事はないのですが...セグメントの開始位置を区別しつつ、cnf と同じフォーマットのコーパスを入力に取れるようにする方法が他に思いつきませんでした。）、せっかく IOB を使うのだから、従来使えていた cnf と同じ素性も使えるようにしてみました。
IOB タグを使い、cnf で使っていたトークン素性を使えるようにするという事は、あるトークン x がセグメントの開始 B になりやすいだとか、また逆に B にはなりにくいという様な情報を学習できます。
さらに、 cnf では素性テンプレートを使い、カレントの位置 i の n 個前のトークンが何であったか等も素性として扱えるため、セグメントの外にあるトークンの情報も扱えるという利点があります。
まさにこれと同じ事をやっているっぽい論文
A Hybrid Markov/Semi-Markov Conditional Random Field for Sequence Segmentation
があったので、名前をお借りして a hybrid markov/semi-markov cnf と名付けてみました。
上の論文自体はまだちゃんと読んで理解したわけではないので、自信を持って同じとも違うとも言えないのですが、ちらっと読んだところ、私と同じくセグメント素性にトークン素性を追加したのではなく、セグメント素性をトークン素性の和で表しているように見えます。
そうするとコストの計算時に、トークン素性のコストを予め計算しておいて、各セグメントで共通する部分についてはキャッシュを使うという事が出来るようになります。
ちゃんと読んでからでないと断言はできませんが、あってたらそこが私の作った学習器との違いになりそうです。
とりあえず、コードは以前 cnf をアップした所にそのまま追加しました。
http://code.google.com/p/cnf/
semi* が追加部分です。
素性テンプレートは以下のようになっています。
トークン素性については以前と変わりありません。
の部分が追加したテンプレートの仕様で、
S:%s[ length , gate function に与える重みの初期値 ]
が unigram segment のテンプレート、
T:%s[ length , gate function に与える重みの初期値 ]
が bigram segment のテンプレートになります。
length は考慮するセグメントの最大長です。
S はセグメントの S ですが、 T は良い頭文字が思いつかなかったので、とりあえず T にしました。
トークン素性の bigram の B 同様、セグメントのラベル遷移だけを素性にする場合は T だけ書いておけば OK です。
ソースを落としてきたら、
みたいにして使います。
コーパスは従来の系列ラベリングで使えた物は大抵使えるんじゃないかと思いますが、セグメンテーションに適した手法なのでどうせならちゃんとしたタスクのコーパスで評価したいところですね。
調べてみた所、公開されている semi-crf の実装もあまり無いようですし、素性テンプレートの使える semi-markov cnf はそこそこ貴重なのではないかと自分で思ってみたりしています。
今回から、学習器の学習率のスケジューリングとL1 正則化は
Stochastic Gradient Descent Training for L1-regularized Log-linear Models with Cumulative Penalty
の手法に変更しました。
@unnonouno さんのブログ
http://unnonouno.blogspot.com/2010/02/l1-sgd.html
で紹介されていた手法で、説明が秀逸ですw
[2010.7.21 追記]
template の仕様を変更しました。
以前の仕様では、unigram segment と bigram segment で各1つずつしかテンプレートを記述できなかったのですが、複数のテンプレートを記述できるようにしました。
それに伴って、
S:%s[ length , gate function に与える重みの初期値 ]
から、
S:%s[ length, col, 重みの初期値]
と、col を指定できるようにテンプレートのマクロも変更しています。
おまけで、新しい template の例を以下に。
semi-crf は入力系列に対する最適なセグメンテーションを学習する学習器です。
crf が P(Y|X, W) を求めていたのに対して、semi-crf では P(S|X, W) を推定します。
文章で書くと、2つの違いは
* crf では入力系列 X に対して生成されうるすべてのラベル系列 Y' に対して、正解系列とそれ以外の系列を弁別するよう学習する
* semi-crf では入力系列 X に対して生成されうるすべてのセグメンテーション結果 S' に対して、正解のセグメンテーションとそれ以外のセグメンテーションを弁別するよう学習する
ということになります。
semi-crf では入力系列が与えられた際に、t_j , ... , u_j までを1つのセグメント s_j として、1つのラベルを割り当てます。
この時、t_j 、 u_j は常に以下を満たす必要があります。
1 <= t_j <= u_j <= |S| , t_{j+1} = u_j + 1
これは、セグメントの終了位置と次のセグメントの開始位置は一致するよということです。
で、crf で使われていた素性関数 f の代わりに segment feature function g ( g^k(y_j,y_{j-1},X,t_j,u_j) )を使います。
といっても、あまり大きく変わっているという事は無くて、crf の素性関数がトークン x_i とラベル y_i のペアに対して 1 or 0 を返していたのが、セグメント s_i とラベル y_i のペアに変わった程度です。
大きな違いは以上なのですが、パラメータ推定をする際のラティスの形が違います。
crf ではラティスはきれいな格子でしたが、semi-crf では各位置 i について、i で始まるセグメントと i で終了するセグメントのみ接続を許すため、結構複雑な形をしています。
ただ、これは形態素解析の物と同じなので 自然言語処理-基礎と応用- の1章、形態素解析の 1.1.5 節あたりを読むと分かりやすいと思います。
パラメータ推定は crf と同じ方法が使えます。
とりあえずここまで理解したところで、cnf で作ってみよう！ と考えて作り始めたのですが、いざ作ってみて失敗に気づきました。
semi-crf ではセグメント単位で同じラベルをふるため、ラベル付けは IOB とかじゃ無く、トークン列に対して AAA、BBB のようにつければ良いやと思っていたのですが、いざ動かす所で、正解の与え方ではまりました。
というのも、AA, AA, AA というラベル列を正解コーパスでつけたつもりでも、AAA, AAA、AAAAAA のどれとも区別が出来なくなってました。。
それで、方向を修正し、結局ラベル付けは IOB タグに限定する事にして（cnf ではそんな事はないのですが...セグメントの開始位置を区別しつつ、cnf と同じフォーマットのコーパスを入力に取れるようにする方法が他に思いつきませんでした。）、せっかく IOB を使うのだから、従来使えていた cnf と同じ素性も使えるようにしてみました。
IOB タグを使い、cnf で使っていたトークン素性を使えるようにするという事は、あるトークン x がセグメントの開始 B になりやすいだとか、また逆に B にはなりにくいという様な情報を学習できます。
さらに、 cnf では素性テンプレートを使い、カレントの位置 i の n 個前のトークンが何であったか等も素性として扱えるため、セグメントの外にあるトークンの情報も扱えるという利点があります。
まさにこれと同じ事をやっているっぽい論文
A Hybrid Markov/Semi-Markov Conditional Random Field for Sequence Segmentation
があったので、名前をお借りして a hybrid markov/semi-markov cnf と名付けてみました。
上の論文自体はまだちゃんと読んで理解したわけではないので、自信を持って同じとも違うとも言えないのですが、ちらっと読んだところ、私と同じくセグメント素性にトークン素性を追加したのではなく、セグメント素性をトークン素性の和で表しているように見えます。
そうするとコストの計算時に、トークン素性のコストを予め計算しておいて、各セグメントで共通する部分についてはキャッシュを使うという事が出来るようになります。
ちゃんと読んでからでないと断言はできませんが、あってたらそこが私の作った学習器との違いになりそうです。
とりあえず、コードは以前 cnf をアップした所にそのまま追加しました。
http://code.google.com/p/cnf/
semi* が追加部分です。
素性テンプレートは以下のようになっています。
# Unigram Segment
S:%s[15,0]
# Bigram Segment
T
# Unigram
U00:%x[-2,0,0]
U01:%x[-1,0,0]
U02:%x[0,0,0]
U03:%x[1,0,0]
U04:%x[2,0,0]
U05:%x[-1,0,0]/%x[0,0,0]
U06:%x[0,0,0]/%x[1,0,0]
U10:%x[-2,1,0]
U11:%x[-1,1,0]
U12:%x[0,1,0]
U13:%x[1,1,0]
U14:%x[2,1,0]
U15:%x[-2,1,0]/%x[-1,1,0]
U16:%x[-1,1,0]/%x[0,1,0]
U17:%x[0,1,0]/%x[1,1,0]
U18:%x[1,1,0]/%x[2,1,0]
U20:%x[-2,1,0]/%x[-1,1,0]/%x[0,1,0]
U21:%x[-1,1,0]/%x[0,1,0]/%x[1,1,0]
U22:%x[0,1,0]/%x[1,1,0]/%x[2,1,0]
# Bigram
B
トークン素性については以前と変わりありません。
# Unigram Segment
S:%s[15,0]
# Bigram Segment
T
の部分が追加したテンプレートの仕様で、
S:%s[ length , gate function に与える重みの初期値 ]
が unigram segment のテンプレート、
T:%s[ length , gate function に与える重みの初期値 ]
が bigram segment のテンプレートになります。
length は考慮するセグメントの最大長です。
S はセグメントの S ですが、 T は良い頭文字が思いつかなかったので、とりあえず T にしました。
トークン素性の bigram の B 同様、セグメントのラベル遷移だけを素性にする場合は T だけ書いておけば OK です。
ソースを落としてきたら、
# make
# ./src/semicnflearn src/semitemplate data/conll2000/train.txt semicnf.save
# ./src/semicnftagger src/semitemplate semicnf.save data/conll2000/test.txt
みたいにして使います。
コーパスは従来の系列ラベリングで使えた物は大抵使えるんじゃないかと思いますが、セグメンテーションに適した手法なのでどうせならちゃんとしたタスクのコーパスで評価したいところですね。
調べてみた所、公開されている semi-crf の実装もあまり無いようですし、素性テンプレートの使える semi-markov cnf はそこそこ貴重なのではないかと自分で思ってみたりしています。
今回から、学習器の学習率のスケジューリングとL1 正則化は
Stochastic Gradient Descent Training for L1-regularized Log-linear Models with Cumulative Penalty
の手法に変更しました。
@unnonouno さんのブログ
http://unnonouno.blogspot.com/2010/02/l1-sgd.html
で紹介されていた手法で、説明が秀逸ですw
[2010.7.21 追記]
template の仕様を変更しました。
以前の仕様では、unigram segment と bigram segment で各1つずつしかテンプレートを記述できなかったのですが、複数のテンプレートを記述できるようにしました。
それに伴って、
S:%s[ length , gate function に与える重みの初期値 ]
から、
S:%s[ length, col, 重みの初期値]
と、col を指定できるようにテンプレートのマクロも変更しています。
おまけで、新しい template の例を以下に。
# Unigram Segment
S01:%s[15,0,0]
S02:%s[15,1,0]
# Bigram Segment
T
# Unigram
U00:%x[-2,0,0]
U01:%x[-1,0,0]
U02:%x[0,0,0]
U03:%x[1,0,0]
U04:%x[2,0,0]
U05:%x[-1,0,0]/%x[0,0,0]
U06:%x[0,0,0]/%x[1,0,0]
U10:%x[-2,1,0]
U11:%x[-1,1,0]
U12:%x[0,1,0]
U13:%x[1,1,0]
U14:%x[2,1,0]
U15:%x[-2,1,0]/%x[-1,1,0]
U16:%x[-1,1,0]/%x[0,1,0]
U17:%x[0,1,0]/%x[1,1,0]
U18:%x[1,1,0]/%x[2,1,0]
U20:%x[-2,1,0]/%x[-1,1,0]/%x[0,1,0]
U21:%x[-1,1,0]/%x[0,1,0]/%x[1,1,0]
U22:%x[0,1,0]/%x[1,1,0]/%x[2,1,0]
# Bigram
B
2010年6月13日日曜日
グラフラプラシアンで推薦
以前縁あって小町さんと一緒に仕事をさせてもらい論文に名前を載せてもらったのですが、会社だけでなく自宅でもちょっと使いたいなーということもあり、実装してみることにしました。
参考にしたのは以下の論文です。
ラプラシアンラベル伝播による検索クリックスルーログからの意味カテゴリ獲得
元論文と違うのは、インスタンス-パターン行列の要素を単純な頻度から別の尺度に変えている点です。
元々そのまんま実装してみたところ、非常にレアな場合なのですが、ジェネリックパターン1つのみと共起するようなインスタンスがあった場合に、これが上位に出やすくなるという問題が発生し、どうにかできないかなと模索していたところ、小町さんからアドバイスを頂き、それを基に手を加えています。
とりあえず動作検証のためにMovieLens Data Setsを使って実験してみました。
最初にデータのフォーマットをツールの入力形式へ変更。
perl swap.pl < ../ml-data/u.data | LC_ALL='C' sort > mlens.data
入力ファイルのフォーマットは次のようにTSV形式で＜インスタンス-パターン-評価値＞として、インスタンスでソートしておきます。
で、実行。
./quetchup -f mlens.data -o test.save
デフォルトでの実行は alpha = 0.0001, iteration = 10, エッジのカットに使う閾値 0.1 としています。
./quetchup -c sample.txt -r test.save | head -50
sample.txt
白雪姫がお気に入りとすると、以下のようにディズニー系の映画を上位にレコメンドすることができます。
ですが、よく見るとスタートレックとかが入り込んでいます。
そこで、次のように嫌いなものリストを用意します。
sample2.txt
./quetchup -c sample.txt -n sample2.txt -r test.save | head -50
と、無事スタートレックとか消えてくれました（スタートレック好きな人ごめんなさい！）。
ちなみに、head で出力を絞っているので見えていませんが、出力にはスコア付きのパターンリストも含まれます。
これで、どんなパターンで各インスタンスがつながっているのかも分かりやすいです。
お気に入りリストが1つしかなくても十分な精度が出そうだし、複数あっても無問題。
嫌いなものあればそれに近いものを推薦候補から外しやすい。
この手法で推薦とか悪くないなと思う今日この頃な訳です。
ソースコードを公開したいのですが、Quetchup アルゴリズムは一応特許が取られているはずなので、それがどんな内容か分からないと不安だったりします。。
参考にしたのは以下の論文です。
ラプラシアンラベル伝播による検索クリックスルーログからの意味カテゴリ獲得
元論文と違うのは、インスタンス-パターン行列の要素を単純な頻度から別の尺度に変えている点です。
元々そのまんま実装してみたところ、非常にレアな場合なのですが、ジェネリックパターン1つのみと共起するようなインスタンスがあった場合に、これが上位に出やすくなるという問題が発生し、どうにかできないかなと模索していたところ、小町さんからアドバイスを頂き、それを基に手を加えています。
とりあえず動作検証のためにMovieLens Data Setsを使って実験してみました。
最初にデータのフォーマットをツールの入力形式へ変更。
perl swap.pl < ../ml-data/u.data | LC_ALL='C' sort > mlens.data
#!/opt/local/bin/perl
# user id | item id | rating | timestamp
use warnings;
use strict;
my %dic;
my $file = "../ml-data/u.item";
open (F, $file) or die;
while ()
{
chomp;
my @t = split(/\|/);
unless (defined $dic{$t[0]})
{
$dic{$t[0]} = $t[1];
}
}
close (F);
while (<>)
{
chomp;
my @tokens = split(/\t/);
print $dic{$tokens[1]}, "\t$tokens[0]\t$tokens[2]\n";
}
入力ファイルのフォーマットは次のようにTSV形式で＜インスタンス-パターン-評価値＞として、インスタンスでソートしておきます。
'Til There Was You (1997) 152 4
'Til There Was You (1997) 178 3
'Til There Was You (1997) 223 1
'Til There Was You (1997) 299 2
'Til There Was You (1997) 342 1
'Til There Was You (1997) 416 3
'Til There Was You (1997) 530 2
'Til There Was You (1997) 532 3
'Til There Was You (1997) 782 2
...
で、実行。
./quetchup -f mlens.data -o test.save
デフォルトでの実行は alpha = 0.0001, iteration = 10, エッジのカットに使う閾値 0.1 としています。
./quetchup -c sample.txt -r test.save | head -50
sample.txt
Snow White and the Seven Dwarfs (1937)
白雪姫がお気に入りとすると、以下のようにディズニー系の映画を上位にレコメンドすることができます。
[Instance]
Snow White and the Seven Dwarfs (1937):0.999801
Cinderella (1950):5.19753e-07
Fantasia (1940):5.06323e-07
Dumbo (1941):4.17243e-07
Beauty and the Beast (1991):3.94127e-07
Aladdin (1992):3.88318e-07
Mary Poppins (1964):3.71035e-07
Jungle Book, The (1994):3.67852e-07
Sword in the Stone, The (1963):3.67118e-07
Lion King, The (1994):3.64219e-07
Pinocchio (1940):3.60772e-07
Robin Hood: Prince of Thieves (1991):3.5665e-07
Aristocats, The (1970):3.51889e-07
Fox and the Hound, The (1981):3.34974e-07
Star Trek VI: The Undiscovered Country (1991):3.24104e-07
Alice in Wonderland (1951):3.19468e-07
Homeward Bound: The Incredible Journey (1993):2.97177e-07
Sound of Music, The (1965):2.97109e-07
Walk in the Sun, A (1945):2.96338e-07
Very Natural Thing, A (1974):2.96338e-07
Bedknobs and Broomsticks (1971):2.93205e-07
Star Trek IV: The Voyage Home (1986):2.90837e-07
Gone with the Wind (1939):2.89412e-07
Nightmare Before Christmas, The (1993):2.85079e-07
Star Trek: Generations (1994):2.8431e-07
Winnie the Pooh and the Blustery Day (1968):2.82242e-07
101 Dalmatians (1996):2.81954e-07
Swiss Family Robinson (1960):2.6257e-07
Aladdin and the King of Thieves (1996):2.54459e-07
Grease (1978):2.53118e-07
African Queen, The (1951):2.52688e-07
Ghost (1990):2.52399e-07
Hunt for Red October, The (1990):2.4501e-07
Miracle on 34th Street (1994):2.42372e-07
Home Alone (1990):2.4004e-07
Top Gun (1986):2.35881e-07
Abyss, The (1989):2.35302e-07
Little Women (1994):2.33329e-07
Old Yeller (1957):2.32874e-07
Pocahontas (1995):2.31783e-07
Mrs. Doubtfire (1993):2.31762e-07
Jumanji (1995):2.31231e-07
Pete's Dragon (1977):2.29198e-07
Christmas Carol, A (1938):2.28481e-07
New York Cop (1996):2.28269e-07
Maltese Falcon, The (1941):2.26072e-07
Goofy Movie, A (1995):2.25177e-07
James and the Giant Peach (1996):2.23962e-07
Star Trek: The Motion Picture (1979):2.23878e-07
ですが、よく見るとスタートレックとかが入り込んでいます。
そこで、次のように嫌いなものリストを用意します。
sample2.txt
Walk in the Sun, A (1945)
Star Trek: Generations (1994)
./quetchup -c sample.txt -n sample2.txt -r test.save | head -50
[Instance]
Snow White and the Seven Dwarfs (1937):0.9998
Cinderella (1950):2.96839e-07
Fantasia (1940):2.81548e-07
Dumbo (1941):2.54658e-07
Aladdin (1992):2.4008e-07
Lion King, The (1994):2.30137e-07
Beauty and the Beast (1991):2.10186e-07
I Don't Want to Talk About It (De eso no se habla) (1993):2.01313e-07
Jurassic Park (1993):1.96396e-07
Coldblooded (1995):1.95535e-07
Robin Hood: Prince of Thieves (1991):1.74025e-07
Fox and the Hound, The (1981):1.71627e-07
African Queen, The (1951):1.69855e-07
Sound of Music, The (1965):1.69781e-07
Mary Poppins (1964):1.64425e-07
Wizard of Oz, The (1939):1.6242e-07
Swan Princess, The (1994):1.56013e-07
Homeward Bound: The Incredible Journey (1993):1.53814e-07
Mamma Roma (1962):1.52174e-07
Good Man in Africa, A (1994):1.50378e-07
Some Like It Hot (1959):1.46221e-07
Citizen Kane (1941):1.42743e-07
Anne Frank Remembered (1995):1.40035e-07
Dances with Wolves (1990):1.36046e-07
It's a Wonderful Life (1946):1.3514e-07
Hunt for Red October, The (1990):1.3501e-07
Young Frankenstein (1974):1.34608e-07
Rebecca (1940):1.32806e-07
Winnie the Pooh and the Blustery Day (1968):1.30989e-07
It Happened One Night (1934):1.30592e-07
M (1931):1.30189e-07
Maltese Falcon, The (1941):1.29532e-07
Lassie (1994):1.28929e-07
Foreign Correspondent (1940):1.28287e-07
Grease (1978):1.2729e-07
Alice in Wonderland (1951):1.23873e-07
Quiz Show (1994):1.22688e-07
Singin' in the Rain (1952):1.22222e-07
Sword in the Stone, The (1963):1.18912e-07
Lawrence of Arabia (1962):1.15236e-07
To Kill a Mockingbird (1962):1.14732e-07
Swiss Family Robinson (1960):1.14698e-07
Drop Dead Fred (1991):1.13925e-07
Vertigo (1958):1.13817e-07
Top Gun (1986):1.12622e-07
Bonnie and Clyde (1967):1.12262e-07
Blue Angel, The (Blaue Engel, Der) (1930):1.10351e-07
Annie Hall (1977):1.09473e-07
Turning, The (1992):1.07052e-07
と、無事スタートレックとか消えてくれました（スタートレック好きな人ごめんなさい！）。
ちなみに、head で出力を絞っているので見えていませんが、出力にはスコア付きのパターンリストも含まれます。
これで、どんなパターンで各インスタンスがつながっているのかも分かりやすいです。
お気に入りリストが1つしかなくても十分な精度が出そうだし、複数あっても無問題。
嫌いなものあればそれに近いものを推薦候補から外しやすい。
この手法で推薦とか悪くないなと思う今日この頃な訳です。
ソースコードを公開したいのですが、Quetchup アルゴリズムは一応特許が取られているはずなので、それがどんな内容か分からないと不安だったりします。。
2010年3月29日月曜日
gibbs sampling
3月に入って netwalker を購入しました。
電車通勤の時間を利用してコードを書きたいなと思っていたのですが、なかなかよろしい感じです。
といってもキーボードは使いやすいとは言い難く、長いコードは書きたくないですが...
ということで、netwalker で作成したプログラム第一号を公開します。
第一号は gibbs sampling を使ったモチーフ抽出アルゴリズムの実装です。
前々から gibbs sampling について調べてはいたのですが、適度な練習問題が無く実装はしていませんでした。
今月になって購入したバイオインフォマティクスの数理とアルゴリズムにちょうど良い例があったので、それを実装してみました。
そもそもモチーフって何という話ですが、一言で言うとタンパク質の配列パターンの事を指すそうです。
私は専門ではないので詳しい事は分かりませんが、ここでは与えられた複数のタンパク質配列から、特徴的な配列パターンを見つけてくる事をモチーフ抽出と言ってよさそうです。
具体的にどうモチーフを抽出するかですが、載っている例では、各タンパク質A、C、T、Gの出現確率をP(A)、P(C)、P(T)、P(G)、部分配列中の位置 j に各タンパク質が現れる確率を f_j(A)、f_j(C)、f_j(T)、f_j(G)とした時に、相対エントロピーを最大化するように部分配列を選んでくる事でモチーフを抽出します。
これをローカルマルチプルアラインメントと呼ぶそうですが、この問題はまともにやると NP 困難なので、gibbs sampling を使って効率的に解くよ、というお話しになっています。
という事で、アルゴリズムは以下の通り。
1. 各入力配列 s_k から, 長さ L の部分配列 t_k を一様ランダムに選択 (k = 1,...,n)
2. 入力配列から, s_i を一様ランダムに選択
3. f'_j(a) を t_i 意外の t_k から計算したモチーフ領域の j 番目の列に置ける文字 a の頻度とする ( f'_j(a) = (|{k|t_k[j] = a, k != i}|) / (n-1))
4. s_i から長さ L の部分配列 t'_i を, Πj f'_j(t'_i[j])/p(t'_i[j]) に比例する確率でランダムに選択
5. t_i を t'_i で置き換える
6. 2-5 を十分な回数繰り返す
ステップ 3. の頻度は確率の間違えかな？と思うのですが、一応そのまま書いておきます。
ステップ 4. では求めた事後分布に従って新しくモチーフ候補を選択し直すわけですが、ここはランダムで選択しないと初期値から動かなくなるので、まず計算したスコアの配列を確率分布として使えるように総和が1となるように正規化してから、PRMLの下巻に紹介されている棄却サンプリングを使いました。
提案分布には一様分布を使っています。
コードと実行例は以下。
サンプルデータ1
サンプルデータ2
実行例1
実行例2
上手い事特徴的な文字列を抽出できる事が確認できました。
電車通勤の時間を利用してコードを書きたいなと思っていたのですが、なかなかよろしい感じです。
といってもキーボードは使いやすいとは言い難く、長いコードは書きたくないですが...
ということで、netwalker で作成したプログラム第一号を公開します。
第一号は gibbs sampling を使ったモチーフ抽出アルゴリズムの実装です。
前々から gibbs sampling について調べてはいたのですが、適度な練習問題が無く実装はしていませんでした。
今月になって購入したバイオインフォマティクスの数理とアルゴリズムにちょうど良い例があったので、それを実装してみました。
そもそもモチーフって何という話ですが、一言で言うとタンパク質の配列パターンの事を指すそうです。
私は専門ではないので詳しい事は分かりませんが、ここでは与えられた複数のタンパク質配列から、特徴的な配列パターンを見つけてくる事をモチーフ抽出と言ってよさそうです。
具体的にどうモチーフを抽出するかですが、載っている例では、各タンパク質A、C、T、Gの出現確率をP(A)、P(C)、P(T)、P(G)、部分配列中の位置 j に各タンパク質が現れる確率を f_j(A)、f_j(C)、f_j(T)、f_j(G)とした時に、相対エントロピーを最大化するように部分配列を選んでくる事でモチーフを抽出します。
これをローカルマルチプルアラインメントと呼ぶそうですが、この問題はまともにやると NP 困難なので、gibbs sampling を使って効率的に解くよ、というお話しになっています。
という事で、アルゴリズムは以下の通り。
1. 各入力配列 s_k から, 長さ L の部分配列 t_k を一様ランダムに選択 (k = 1,...,n)
2. 入力配列から, s_i を一様ランダムに選択
3. f'_j(a) を t_i 意外の t_k から計算したモチーフ領域の j 番目の列に置ける文字 a の頻度とする ( f'_j(a) = (|{k|t_k[j] = a, k != i}|) / (n-1))
4. s_i から長さ L の部分配列 t'_i を, Πj f'_j(t'_i[j])/p(t'_i[j]) に比例する確率でランダムに選択
5. t_i を t'_i で置き換える
6. 2-5 を十分な回数繰り返す
ステップ 3. の頻度は確率の間違えかな？と思うのですが、一応そのまま書いておきます。
ステップ 4. では求めた事後分布に従って新しくモチーフ候補を選択し直すわけですが、ここはランダムで選択しないと初期値から動かなくなるので、まず計算したスコアの配列を確率分布として使えるように総和が1となるように正規化してから、PRMLの下巻に紹介されている棄却サンプリングを使いました。
提案分布には一様分布を使っています。
コードと実行例は以下。
#!/usr/local/bin/perl
use warnings;
use strict;
my $L = 5;
my @dataset;
my %wordset;
my $words = 0;
while(<>)
{
chomp;
my $line = $_;
my @sequence = split(//,$line);
push @dataset,\@sequence;
}
# init words and probs
$words = &setwords(\%wordset, \@dataset);
&initwordprobs($words,\%wordset);
my $datasize = @dataset;
print "dataset\n";
foreach my $ar (@dataset)
{
print "@$ar\n";
}
# init motif set
my @pos;
foreach my $ar (@dataset)
{
my $p = int rand(@$ar-$L);
push @pos, $p;
#print $p,"\n";
}
# output first
print "init\n";
for (my $i = 0; $i < @dataset; $i++)
{
my $ar = $dataset[$i];
my $p = $pos[$i];
for (my $j = $p; $j < $p+$L; $j++)
{
print "$ar->[$j] ";
}
print "\n";
}
# iterration
for (my $iter = 0; $iter < 50; $iter++)
{
#print "epoch: $iter\n";
my $id = int rand($datasize);
my %wordsfreq;
# wordset * position
# +---+-------------+
# | |a|b|c|d|e|f|g|
# +---+-------------+
# | 0 | |
# | 1 | |
# | 2 | frequency |
# | . | |
# | L | |
# +---+-------------+
&initwordsfreq(\%wordset,$L,\%wordsfreq);
# count word frequency
&countwordfreq($id,\@dataset,\%wordsfreq,\@pos);
# pick up new motif i
my $ar = $dataset[$id];
$pos[$id] = &pickup($id,$ar,$L,\%wordsfreq, \%wordset);
}
# output
print "result\n";
for (my $i = 0; $i < @dataset; $i++)
{
my $ar = $dataset[$i];
my $p = $pos[$i];
for (my $j = $p; $j < $p+$L; $j++)
{
print "$ar->[$j] ";
}
print "\n";
}
# subroutine
sub pickup
{
my $id = shift;
my $ar = shift;
my $L = shift;
my $wordsfreqhr = shift;
my $wordsethr = shift;
my @scores;
my $total = 0;
for (my $p = 0; $p <= @$ar-$L; $p++)
{
my $score = 1;
my $index = 0;
for (my $j = $p; $j < $p+$L; $j++)
{
my $w = $ar->[$j];
$score *= $wordsfreqhr->{$index}->{$w}/$wordsethr->{$w};
$index++;
}
push @scores, $score;
$total += $score;
}
return int rand(@scores) if $total == 0;
for(my $i = 0; $i < @scores; $i++)
{
$scores[$i] /= $total;
}
#print "@scores\n";
my $flg = 1;
while($flg)
{
my $npos = int rand(@scores);
my $s = rand(1);
if ($s <= $scores[$npos])
{
$flg = 0;
return $npos;
}
}
}
sub countwordfreq
{
my $id = shift;
my $datasetar = shift;
my $wordsfreqhr = shift;
my $posar = shift;
for (my $i = 0; $i < @$datasetar; $i++)
{
next if ($i == $id);
my $ar = $datasetar->[$i];
for (my $j = $posar->[$i]; $j < $L+$posar->[$i]; $j++)
{
my $p = $j-$posar->[$i];
my $w = $ar->[$j];
$wordsfreqhr->{$p}->{$w}++;
}
}
for (my $l = 0; $l < $L; $l++)
{
foreach my $k (keys %{$wordsfreqhr->{$l}})
{
$wordsfreqhr->{$l}->{$k} /= (@$datasetar-1);
}
}
}
sub initwordsfreq
{
my $wordsethr = shift;
my $L = shift;
my $wordfhr = shift;
for (my $i = 0; $i < $L; $i++)
{
my %freqonpos;
foreach my $k (keys %$wordsethr)
{
$freqonpos{$k} = 0;
}
$wordfhr->{$i} = \%freqonpos;
}
}
sub initwordprobs
{
my $words = shift;
my $wordsethr = shift;
foreach my $k (keys %$wordsethr)
{
$wordsethr->{$k} = 1/$words;
}
}
sub setwords
{
my $wordshr = shift;
my $dataar = shift;
my $words = 0;
foreach my $ar (@$dataar)
{
foreach my $w (@$ar)
{
unless (defined $wordshr->{$w})
{
$wordshr->{$w} = defined;
$words++;
}
}
}
return $words;
}
サンプルデータ1
c a t c g a a a a a a
a a a a c a t c g a a
a a c a t c g a a a a
a a a a a a c a t c g
サンプルデータ2
t c a g a g a t t c g t a g
t t c a t t c g g g c g t
c g a t t c g c g a c t c
t g a a t t c g g a a
実行例1
# perl gibbs.pl < sample1.txt
dataset
c a t c g a a a a a a
a a a a c a t c g a a
a a c a t c g a a a a
a a a a a a c a t c g
init
g a a a a
a a a c a
a t c g a
a a a a a
result
c a t c g
c a t c g
c a t c g
c a t c g
実行例2
# perl gibbs.pl < sample2.txt
dataset
t c a g a g a t t c g t a g
t t c a t t c g g g c g t
c g a t t c g c g a c t c
t g a a t t c g g a a
init
t c a g a
c a t t c
t c g c g
a a t t c
result
a t t c g
a t t c g
a t t c g
a t t c g
上手い事特徴的な文字列を抽出できる事が確認できました。
登録:
投稿 (Atom)
投稿
