Probabilistic Graphical Models(pgm) - 1466006216

1：karino2：2016/06/16(木) 00:56:56: Courseraのco-founderのコース。

ベイジアンネットワークのコース。
ベイジアンネットワークは各ノードがベイズの推論となっていて、エッジが条件付き確率の条件を表すような物。

最初このコースを見た時はなんのコースかすら分からなかったですが、MCMC勉強した後にみると凄くMCMCと関連したコースである事に気づく。
8：karino2：2016/06/25(土) 00:56:04: ここらでクイズを解くか、という事でWeek1のクイズに挑む。
Q 6まではすらすらと進むが、Q7でちょっと手が止まる。良いね。
9：karino2：2016/06/25(土) 01:04:30: Q8のI-mapで分からなくなったので、スライド見直し。
Pがfactorize over G: Pの同時確率が、各Xiのグラフ上の親の元での条件付き確率の積で表される事。
10：karino2：2016/06/25(土) 02:20:58: やはりI-mapが分かっていない。
動画を見直す。

Gから考えられる独立のノードの組み全てが、Pの上でも独立ならGはPのI-map、という事か。
別にPはGで表されるベイジアンネットワークとは言っていない。
どういう違いがありうるのか？良く分からないな。
11：karino2：2016/06/25(土) 02:30:34: 動画2-5を見直した。なるほど。
Pがグラフから考えられるD-separationからくる独立をすべて満たしていれば、GはI-mapだ、と言えると。
この時Pの独立な組み合わせ全てをGが表現している必要は無い。
また、独立以外の要素については何も言及していない。
つまりPがGのベイジアンネットワークである事は要請していない。
12：karino2：2016/06/26(日) 23:00:12: I-mapは大分理解したので動画6-5を見ていた。
I-equivalentなるものが出てきて、なるほどぉ、とか思うなど。
13：karino2：2016/06/29(水) 00:33:59: 6-6と6-7はいまいち何が言いたいのか良く分からないな。
6-6はLog linearモデル。これが一般的なモデルなのはまぁいいんだが、PGMとの関係が分からん。
6-7のLog linearモデルのshared featureも、そういう事はあるのだろうがだから何？という印象。

これ、なんかこう、こういうモデルを一般的に解くアプリとかあって設定するだけで使えたりするのかなぁ。
それだと意味は分かるが。

そうでなければ多分何も理解出来てない。
6-7のsummaryのテンプレートを使い回せるとかパメータと構造をMN内、及び異なるMN間で使い回せる、というのもいまいち分かってないので、たぶん理解出来てない気がするなぁ。

もうしばらく先に進めてみて、分かって無さそうだったら元に戻ろう。
14：karino2：2016/07/17(日) 08:46:50: 適当に寝る前とかに動画を見続けて、現在9.1のBelief Propagationまで来た。
この辺はやれば出来そうな気がする。
Variable Eliminationとかこの辺よりも、もっと基礎的な所の動画を見直した方が良いかもなぁ。
とりあえず分かる所まで見続ける予定だが。
15：karino2：2016/07/18(月) 16:36:45: クラスタグラフのメモ

- Family Preservation
- Running Interception Property
16：karino2：2016/07/20(水) 21:00:37: Loopy BPまで見終わり。
なんかこの話は可能性を感じさせていいね。
17：名無しさん：2016/12/08(木) 04:25:01: 12/5から再開したようなので、いい機会と再履修。

Factorは、引数のランダム変数が特定の値を取った時に、そこから実数値への関数、という事でいいのかな。
18：karino2：2016/12/08(木) 05:16:16: url貼っとく。

https://www.coursera.org/specializations/probabilistic-graphical-models
19：karino2：2016/12/12(月) 00:23:51: 今の所お金払って受ける気は無いが、プログラミングエクササイズは途中までは自習するかな。
現在Week 1の動画を見終わってプログラム課題をやっている。
問題のストーリー設定が無駄に細かくて良いね。
20：karino2：2016/12/12(月) 01:51:08: SamIamで課題のネットワークを作成。面白いね、これ。
21：名無しさん：2016/12/16(金) 09:51:59: Week 1のプログラム課題を見始めたが、難しいねぇ。
開幕のFactorProductからして何を作りたいか理解するのに一時間くらい時間かかった気がする（Ocvave忘れてたというのもあるが、、、）

ようするに同時確率のテーブルとその演算を行うライブラリを構築しようとしているんだな。
これが最初の課題かよ、、、
22：karino2：2016/12/16(金) 09:54:56: varは順番だけで、実際に値はcardinalityを元に、例えば二次元の時は

1, 1
2, 1
1, 2
2, 2

と自動で振られるのがポイントだな。
23：karino2：2016/12/16(金) 10:08:45: 最初のTutorialの例はこういう事だな。

https://imgur.com/pdxnKyE
24：karino2：2016/12/16(金) 10:17:58: テスト

https://imgur.com/pdxnKyE.jpg
25：karino2：2016/12/17(土) 06:39:55: Factorizeは問題を理解するのに1時間くらいかかったが、理解すればすぐ解ける。
で、次のMarginalizeに挑んでいる。

これは指定されたベクトルの各要素の1:cardな配列を作って、それの一致する要素だけでsumをとればいいのだが、
これをどうOctaveで表現するのかが分からない。
たぶんそのまんまコードで表現出来ると思うのだが…
26：karino2：2016/12/17(土) 06:49:35: そうか、Bのassignから始めて、margiinalizeする部分を1:cardとすれば良いのか。
27：karino2：2016/12/17(土) 07:17:34: 4, (1:3), 5] みたいな感じで[4, 1,5], [4, 2, 5], [4, 3, 5]が作りたい。
やり方は分からないが、[4, 5]を繰り返して、そこに1:3を追加して、A(:, [2, 1, 3])とかすればいいんじゃないか。
28：karino2：2016/12/17(土) 07:25:09: >> v = 1:4
v =

1 2 3 4

>> repelems(v, [1, 2, 3, 4; 2, 2, 2, 2])
ans =

1 1 2 2 3 3 4 4

あと一歩だなぁ。
29：karino2：2016/12/17(土) 07:41:24: >> ass = IndexToAssignment(1:12, [3, 2, 2])
ass =

1 1 1
2 1 1
3 1 1
1 2 1
2 2 1
3 2 1
1 1 2
2 1 2
3 1 2
1 2 2
2 2 2
3 2 2

で、2についてmarginalizeしたいとする。
>> find(ass(:,2) == 1)
ans =

1
2
3
7
8
9

これで1についてのassignを得る事は出来る。
2についても同様に得られるのでfor文を回せば良いのは分かる。
とりあえずfor文で解いてみよう。
30：karino2：2016/12/17(土) 09:42:58: いろいろ考えて、Bのassignmentに対応したAのindexがあれば良い！という結論になったら、
自分で書くパートの前の所ですでに用意してくれてた…
まぁそういう事を自分で気づける、という必要はあるので、無駄では無かろう。

結果はvectorize出来ると思うのだが、少し考えて分からなかったのでやめておく。
31：karino2：2016/12/17(土) 10:55:27: ObserveEvidenceも解ける。よしよし。
たださすがに疲れたのでここらで一旦休憩。
32：karino2：2016/12/17(土) 18:48:56: 課題一通り解けた！
解けてみると難しい所は無いのだが、問題を理解するのがなかなか大変で時間かかった。
33：karino2：2016/12/18(日) 07:08:50: Week2に入る。
Template Modelの動画とクイズ見終わり。
34：karino2：2016/12/22(木) 14:07:25: week2のプログラム課題開始。
最初の小話は要るのか？^^;
35：karino2：2016/12/22(木) 15:08:35: 最初の問題にたどり着いたが、その前にAppendixで遺伝学の勉強だす。
36：karino2：2016/12/22(木) 15:21:07: メンデル的な遺伝のケースを読み終わる。
英単語の勉強としても教養的な背景知識の勉強としても短くて良く書けているな。

自分でコンピュータの本とか幾ら読んでいてもこういう語彙を学ぶ機会は無いよなぁ。
Courseraは留学っぽい効果があるよなぁ、と思う所である。
37：karino2：2016/12/23(金) 03:08:37: non mendel な方まで終わり。
2の時-1で1の時1になるような式をでっちあげてベクトル化に成功。
38：karino2：2016/12/23(金) 08:40:40: alleleFreqの方もベクトル化まで含めて終了。
順調順調♪
39：karino2：2016/12/23(金) 08:43:01: 一つ一つの問題が、仕様を理解するのに相当程度時間がかかる。
そこまで行け実装自体は10行以下なのだが。
40：karino2：2016/12/23(金) 09:25:06: 親がいる場合のgenotypeの計算だん。
ただベクトル化は諦める。
for が4つネストしてて、中2つくらいはベクトル化しようとしたのだが良く分からなかったので。
41：karino2：2016/12/24(土) 10:48:23: さて、次はconstructGeneticNetwork。
genotypeってFfとfFの区別ないんだっけか。
この辺の試してわかった事が数日経つとすぐ忘れるなぁ。
どうにか残しておいた方が良いかもしれん。
42：karino2：2016/12/24(土) 11:07:28: これまで作った関数を呼び出すだけでだいたい出来た。
ベクトル化はできそうだが、将来の自分が読みにくくなるだけっぽいからやめておく。
43：karino2：2016/12/24(土) 11:11:37: クイズを答えるヤツはかったるいので、先にdecoupled BNに進む事にする。
44：karino2：2016/12/24(土) 11:53:04: 愚痴など書く。

http://karino2.livejournal.com/423688.html
45：karino2：2016/12/24(土) 11:58:32: 愚痴という訳でも無いか。
さて、Decoupleを再開しよう。
46：karino2：2016/12/24(土) 20:41:44: constructDecoupledGeneticNetworkの実装が終わった。
なかなか苦労してメモ取りながら進めたのであとで貼る。

もうベクトル化をしようという気すら起こらない、、、
47：karino2：2016/12/25(日) 09:06:16: >> 46

Decoupledを解く。どこから手をつけたもんか。

まずは提供されている物から調べてみよう。
childCopyGivenParentalsFactorを実行してみる。
それを呼ぶ事になるであろうphenotypeGivenCopiesFactor()の呼び出しの引数から、きっと必要な引数は生成出来るに違いないと予想。

シグニチャを見ると以下。

childCopyGivenParentalsFactor(numAlleles, geneCopyVarChild, geneCopyVarOne, geneCopyVarTwo)

numAllelesはまぁいい。Varは単なる変数名で良ければ、とりあえず何を入れても動くはずだ。結果のどこに反映されたか分かりやすいように5くらいからの値を入れておこう。5, 6, 7だな。
とりあえず単純なケースとして、FF, Ff, ffのgenotypeのパターンを突っ込んで何が起こるか見る。この場合numAllelesは2か。

>> childCopyGivenParentalsFactor(2, 5, 6, 7)

ans =
scalar structure containing the fields:
var =
5 6 7
card =
2 2 2
val =
1.00000 0.00000 0.50000 0.50000 0.50000 0.50000 0.00000 1.00000

つまりどういう事だろう？
varには5, 6, 7が入るので、条件付き確率のvarとしてそのまま入っているのだろう。
cardはnumAllelesから生成されるのであろうな。
結果のvalは8。これはどういう事だろう？そもそもcardは何故3じゃないのか？
これまでなら、

FF | 父FF, 母FF
Ff | 父FF, 母FF
ff | 父FF, 母FF
FF | 父Ff, 母FF
Ff | 父Ff, 母FF
ff | 父Ff, 母FF
...

という数値がvalに並ぶのだから、27個の値が並ぶはずである。それが8個だけとなっている。
48：karino2：2016/12/25(日) 09:08:33: >>46-47

このcardが何かは分からないが、それをXとYと置こう。さらにPDFから察するに、これの引数はどちらも父親から来ているとの事。
すると求める答えは

X | 父1X, 父2X
Y | 父1X, 父2X
X | 父1Y, 父2X
Y | 父1Y, 父2X
X | 父1X, 父2Y
Y | 父1X, 父2Y
X | 父1Y, 父2Y
Y | 父1Y, 父2Y

の結果が並んでいるはずだ。
なるほど、つまりgenotypeの片方が入るのだろうな。つまりこうか。

F | 父FF
f | 父FF
F | 父fF
f | 父fF
F | 父Ff
f | 父Ff
F | 父ff
f | 父ff

よし、cardは8になった。次に値の解釈だ。
並べて書くと、

1.00000 F | 父FF
0.00000 f | 父FF
0.50000 F | 父fF
0.50000 f | 父fF
0.50000 F | 父Ff
0.50000 f | 父Ff
0.00000 F | 父FF
1.00000 f | 父FF

なるほど、単純にどちらのalleleがgenotypeから選ばれるかの、場合の数からくる確率だな。

次はchildCopyGivenFreqsFactor()。これはpdfから察するに事前分布から想像される親の情報無しでのalleleの分布だな。
シグニチャを見ると以下。

childCopyGivenFreqsFactor(alleleFreqs, geneCopyVar)

変数名は5を入れて、alleleFreqsに、例えば[0.9, 0.1]を入れてみよう。

>> childCopyGivenFreqsFactor([0.9, 0.1], 5)
ans =
scalar structure containing the fields:
var = 5
card = 2
val =
0.90000
0.10000

つまり父がFFとか持つ時の、その片方のFがどういう確率か、という事だな。
Copyというのはその二つを合わせるとFFになるような、片方を表すのだな。だいたい分かった。

では問題のphenotypeGivenCopiesFactor()を実装してみよう。
49：karino2：2016/12/25(日) 10:24:01: spinal musclar atrophy
脊髄性筋萎縮症
50：karino2：2016/12/25(日) 21:08:57: constructSigmoidPhenotypeFactorの結果が合わないので真面目に考える。

geneが2つでそれぞれにalleleが二つあるとすると、それだけで4通り、それが二つの親からくるのだから8通りになる気がするが、
sampleGeneticNetworkを見ると条件付き確率の条件は4つに見える。

PDFのサンプルを読み直して考えてみる。まずAabBのケースの計算がある。
ここでAaAAとかの場合が無い、というのが上記の8通りにならない理由に思う。

Aaはgene 1からのものだな。で、bBがgene 2からだ。
でもAaは二つのCopyからくる（両親）。
えーと、つまりどういう事だ？
つまり条件付き確率は以下の形式って事か？

Pheno | 親1G1 親2G1 親1G2 親2G2

ではgeno1について考える場合、親1と親2のgeno1がこれで分かる訳だ。
お、理解してきた。alleleWeightsはgenoについての重さなのに自分は親で見ていたんだな。

直してみよう。

合わないなぁ。なんかsampleGeneticNetworkのvarを見ていると、親1が先に並んでいるな。こうか？

Pheno | 親1G1 親1G2 親2G1 親2G2

これで試してみよう。

出来た！
ここからベクトル化して一般化するのは簡単ね。
これでWeek 2のPAはおしまい！
51：karino2：2016/12/26(月) 10:27:11: 次はMRF。
とりあえずpairwise、Gibbs分布、CRFは理解した。
以前も見た動画のはずだがいまいち区別に記憶が無いので、以前より理解は上がってる気がする（忘れてるだけかもしれないが）
そしてクイズの答えが知りたい、、、
52：karino2：2016/12/26(月) 13:01:11: MRFの独立の話とBNとMRFの表現出来る独立が違うという話や、minimal I-mapやPerfect I-mapの話が終わる。
この辺は言ってる事一つ一つは分かるんだが、要点を理解しているかは自信無いなぁ。
53：karino2：2016/12/27(火) 09:25:06: week3のPA開始した。
現在ScoreModel を動かすところまではつつがなく進む。

この課題はモデルが本格的で面白いね。
54：karino2：2016/12/27(火) 10:32:20: tripletを考慮に入れてScore呼んだら帰ってこない。
自分の実装部分はそんなに遅く無いと思うんだが。

グローバル変数を使うように変更したがやはり帰ってこない。
sparseに返す方法があるのかなぁ。
フォーラム見るか。
55：karino2：2016/12/27(火) 10:34:23: フォーラムにもなんも載ってない(´・ω・｀)
56：karino2：2016/12/27(火) 13:54:23: 凄い時間がかかる物っぽい。
自分のBaytrailマシンでは30分以上はかかってた。

さらにChooseTopSimilarityFactorも実装して走らせて、現在待ち中。

Week3の課題は全体的に簡単ですね。
しかもモデルがなかなか本格的で面白い。
唯一の問題はあんまり精度上がらないところか。
単語精度で37％だと、まだ使い物にならないなぁ、という印象。

あとでVisualizeWordの結果も見てみたい。
57：karino2：2016/12/27(火) 14:40:44: Vizualizeはなんかフリーズした。
というか全体的にプロットとかはフリーズするんだよなぁ。
58：karino2：2016/12/28(水) 06:21:16: Decision Makingの動画見始め。
思ったより面白い。
この辺の動画、全然見た記憶無いんだが、、、と以前ダウンロードした動画を見直すと、見つからない。
あれ？取り直したのかな？そうは見えないが。

Utilityもfactorに過ぎないのでfactor productをマージナライズすれば良い、ってのは美しいな。
59：karino2：2016/12/29(木) 01:39:01: Week4の動画見終わり、課題のpdfを見ているが、これは、、、
本格的過ぎて怯んでしまうが、今回も理解してしまえば書く所は少ないのだろうなぁ。

でもこんなの毎週の課題で実装しろとか凄いよな。
このペースで実装していけるなら、実務でも大分いろいろ出来そうだが。
60：karino2：2016/12/30(金) 09:55:49: このコースの課題は、すごく良く出来ていて、PGMという物の実際が凄く良く分かる。
むしろこの課題を中心にコースを作る方が良い気がするくらい。
やはり課題ベースの方が勉強はいいんだろうね。
61：karino2：2016/12/30(金) 10:14:30: SimpleCalcExpectedUtilityの実装が終わる。
クイズはかったるいのでとりあえず飛ばす。
もしあとでわからなくなったら戻ってきてやろう。
62：karino2：2016/12/30(金) 11:50:24: OptimizeMEUのヒントが理解出来ないのとサンプルのoptdrの出力が理解出来ないのでフォーラムで質問。
ヒントと出力が間違ってるだけだと思うのだが、、、
63：karino2：2016/12/30(金) 13:51:25: ヒントを理解して実装したら、サンプルの出力の間違いに気づく。

optdrは

2
1 1.0
2 0.0

だね。
なかなか難しいがたぶん理解出来たと思う。
64：karino2：2016/12/30(金) 16:37:50: 次はJointUtility。
pdf ではPA4で求めたかもしれない云々書いてあるが、これは以前のコースの時のだね。修正漏れだ。
以前のWeek4はinference でBPとかMax-Sum Message Passingとかやってる。

とりあえずFactorSomeを実装してみるか。
65：karino2：2016/12/30(金) 17:37:38: ユニットテストの2を実行したら、そもそもSimpleCalcExpectedUtilityの結果が合ってない。
真面目に考え直す。
66：karino2：2016/12/30(金) 20:40:51: 二つ目のテストを走らせたらSimpleCalcExpectedUtilityがそもそも合ってないので、真面目に考え直す。

まず、PDFに書いてあるケースを再現してみる。

>> PO = struct('var', [2,3], 'card', [2, 2], 'val', [0.15, 0.7, 0.1, 0.05]);
>> PrintFactor(PO)
2 3
1 1 0.150000
2 1 0.700000
1 2 0.100000
2 2 0.050000

>> U = struct('var', [2,3], 'card', [2, 2], 'val', [300, 50, -400, -500]);
>> PrintFactor(U)
2 3
1 1 300.000000
2 1 50.000000
1 2 -400.000000
2 2 -500.000000

まずはここから。
単純に要素ごとに積をとって足せば良いはず。

>> sum(PO.val .* U.val)
ans = 15

自分の理解はあってる気がする。

自分の理解では、全部の同時確率をfactorの積で表し、VEでUのparent以外を全部消して、あとはFactorProductで掛けてvalのsumを求めれば良い、というもの。
この手続きで上記のデータで同じ事をやっても、やはり15が出た。

なんかCPDFromFactorがバグってるんじゃないか？
最初X3が以下だが、

>> PrintFactor(X3)
3 1 2
1 1 1 4.000000
2 1 1 4.000000
1 2 1 1.000000
2 2 1 1.000000
1 1 2 1.000000
2 1 2 1.000000
1 2 2 4.000000
2 2 2 4.000000

>> X3 = CPDFromFactor(X3,3);
>> PrintFactor(X3)
1 2 3
1 1 1 0.500000
2 1 1 0.500000
1 2 1 0.500000
2 2 1 0.500000
1 1 2 0.500000
2 1 2 0.500000
1 2 2 0.500000
2 2 2 0.500000

全部0.5になってしまっている。たぶん各var 3のそれぞれをノーマライズしたいんだと思うが。
ちょっとコードを見直そう。

見直した。ようするにparentのassignごとにvar 3の値の確率になるようにノーマライズする。
つまりassignからvar 3のインデックスを抜いた物を固定して、var 3をcardだけ動かした和をとると1になるようにする。
するともともとvar 3の値が1と2の値はどちらも同じなので、これで合ってる。

なるほど。
67：karino2：2016/12/30(金) 20:41:32: cumulativeなarrayfunc的なのが無いかと思ったが無い模様。

http://stackoverflow.com/questions/12144569/is-there-any-other-built-in-cumulative-function-like-cumsum-cumprod-in-matlab

まじかよ。
68：karino2：2016/12/30(金) 21:18:12: 終わった！頑張った！
Week 4の課題もなかなか大変だったが、良く理解は深まった。
良い課題で、やってて楽しい。ただ時間凄いかかるね。ほぼ丸一日かかった。
69：karino2：2016/12/30(金) 21:27:43: 終わったと思ったが、Value of Perfect Informationがなかなか難しい(^^;
もう一息だね。
70：karino2：2016/12/30(金) 22:40:02: 今度こそ出来た～！これはいい問題。
71：karino2：2016/12/31(土) 13:48:45: Week5見終わり。
FinalExamは解かないけど眺めるだけ眺めた。
そしてI-equivalentとかなんだっけ？とスライドを見直すなど。
グラフから読み取れる独立の構造が同じグラフの事かね。I(G_1)=I(G_2)となるG_1とG_2が定義の模様。

これで3つのSpecializationのうち一つ目が終わった。
一つ目はrepresentationという事で構成されてるが、微妙にプログラム課題ではVEとか出てきて以前の構成を引きずっている印象がある。

ただだいたいは分かるし、基本的には収まりは良いし、何よりこれ全部とるの大変なので途中で一旦チェックポイント的なのがあるのは良い気がする。
ただやはりInferenceまではやらないとご利益が薄いので、2つめもやる必要はあるよねぇ、という事で二つ目に進みます。
72：karino2：2017/01/02(月) 12:46:50: 年もあけたし、忘れないうちにPGMの2つ目のコースを開始する。
ちょっとフライングだけどまぁいいでしょう。

二つ目のコースはinference。
まずはSum-Productで目的の確率分布だけを残すという話。
73：karino2：2017/01/02(月) 13:13:54: MAP推定が出てきたので別の本で最尤推定とMAP推定を復習しておく。
ほとんど両者は同じだが、MAP推定は事後分布を使う所が違いか。
74：karino2：2017/01/02(月) 14:35:51: VEのcomplexityの解析。
スコープのcardinalityはr_kと言ってるが、d^r_kじゃないのかしら？
cardinalityという言葉の意味を私が誤解してるのかしら？
75：karino2：2017/01/03(火) 10:18:10: VEの動画、見終わり。
クイズも5くらいまでやる。
この辺は別段分からない事も無いね。
次次��♪
76：karino2：2017/01/03(火) 10:46:45: メッセージパッシング。
クラスターグラフにファクターを割り当ててメッセージを送り合う。
BPとかsepsetがここで登場。
77：karino2：2017/01/06(金) 01:20:19: クラスターグラフの満たすべき2つの性質。
Family PreservationとRunning Intersection Property。
78：karino2：2017/01/08(日) 09:36:19: Blief Propの動画見終わり。
クイズも半分くらいやる。
次は前回良く分からなくなった因縁のクリークツリーのあたり。
79：karino2：2017/01/08(日) 09:50:05: ツリーの定義を確認しようとスライドを見直したが見つからなかったのでwikipediaのリンクを貼っとく。

https://ja.m.wikipedia.org/wiki/%E6%9C%A8_(%E6%95%B0%E5%AD%A6)

連結でループが無いもの、か。
80：karino2：2017/01/08(日) 16:57:44: なんかClique Treeの動画がオフラインダウンロードされてない。
あれー？しなかったかなぁ。
以前ダウンロードしたファイルはあるのでその動画を見る。

今Clique TreeのCorrectnessを見ているが、最後の例、この説明では不十分な気がするので一時停止して考える。

ある変数をマージナライズする時は、そこから先にその変数が存在しないのは、RIPとClique Treeの定義から分かる。
だからマージナライズする時は全てのfactorが掛けられている事が保証されてる。
81：karino2：2017/01/08(日) 17:15:13: computationに移る。
何故デルタ1_2は最初に計算した値が収束した値なのか？

それはクラスター1はクラスター2からのメッセージ以外を受け取らない。
で、デルタ1_2は2からのメッセージを抜いた物の和でマージナライズするから、2からのメッセージが幾ら来てもデルタ1_2は変わらない。

さて、この考察をクラスタ2にも適用していくのだろうけれど、これが成立する条件をもう少し考えたい。
ようするにデルタ2_3を計算する時に1からしかメッセージが来てない、という事が言えればいいんだろうな。
つまり数学的帰納法でなんか証明するんだろう。

まずエッジが一つのノードから始める。これを1と呼ぶ。
で、その隣のノードは1からのメッセージはこの場で確定する。
この隣のノードを2と呼ぼう。
2は複数のノードに繋がってる可能性があるから、その別のノードからのメッセージはやってくる可能性があるはず。

https://imgur.com/EYI9ofn.jpg

これが多分RIPを使って影響無い、といえるのだろうな。
82：karino2：2017/01/08(日) 17:22:31: 次のConvergence のスライドと説明を見てわかった。
別の枝も必ず端点があるはずなんだよな。
で、目指す方向以外からの全端点からのメッセージを待つと、そこまでの全メッセージがその時点で収束している、という事が示せるのだろう。
それはツリー性から示すのだろうな。

1. 端からのメッセージは確定している
2. 端からi番目のメッセージが確定していれば、i+1番目が確定する

たぶん言えそうだな。じゃあいいや（ぉぃ
83：karino2：2017/01/08(日) 17:24:54: ツリーの深さを葉から考えてparentに上がっていくと考えれば多分いけそうだな。
84：karino2：2017/01/08(日) 20:13:15: クイズの図が出てこない。
まだ始まってないからかなぁ。
85：karino2：2017/01/08(日) 21:15:16: VEとClique TreeのMessage Passingの比較が良く分からない。
ようするに計算すると同じ事やってる、となるという空気は伝わってくるのだが。

redunduntって何が？
86：karino2：2017/01/08(日) 23:29:04: Week2の動画、一通り見終わり。
前回見た時も一つ一つの動画が分からないという感じじゃなかったのだが、そもそもなんでそんな話をしてるか、みたいなのが良く分かってなかった。

今回はスペシャリゼーションとしてコースが分割されて、BPはinference としてやってる、というのが意識しやすくなってストーリーを見失う事が減った気がする。

一部わからない所もあるが、全体としては割としっかり理解出来て居るんじゃないか。

プログラム課題はweek3が締め切りとの事だがどうしようかなぁ。やっとこうかなぁ。
87：karino2：2017/01/09(月) 12:23:59: PAのexact inference やるぞ(ﾟДﾟ)ｺﾞﾙｧ！
88：karino2：2017/01/09(月) 22:11:16: さて、とりあえず最初の課題であるComputeInitialPotentials.mまでたどり着いた。
このコースの課題はこの手の、実際の問題設定を理解するまでが結構辛いんだよなぁ。

で、指示通りload("PA4Sample")してInitPotential.INPUTとInitPotential.RESULTを眺める。
nodesってなんだろう？

コメントを見るとツリーの中のクリークを表すcell arrayだそうで。
クリークを表す、というのは、そのfactorを持つクラスタという風に解釈すれば良いのかな。
まずその路線で実装してみるか。

cell arrayの行成分の意味が分からないな。サンプルも全部1だし。
あ、{}で参照するものなんだっけ、そういえば。

nodes{2}を見たら2, 8と書いてあるので、何も考えずにFactorProduct(factorList(2), factorList(8))してみたが、
そもそもvarの値からして違う。まぁそうだよなぁ。FactorProductは変数はfactorの変数で、クリークのvarはfactorのインデックスを表しているのだから。
全く分かってないがmarginalizeしてみたが答え合わない。そういう問題じゃないよな、やっぱ。

でもPDFの2の2.で、ファクターのスコープがクリークに入るように、って書いてあるよなぁ。C_iとはなんぞや。
89：karino2：2017/01/09(月) 22:13:18: お、このサンプルのINPUT、Figure 2のネットワークか。なるほど。

なんか1と3は単なるfacotrListのproductで一致してるなぁ。なんで2は一致しないのか。
それはfacotrList(2)の変数に1が含まれているからじゃないか。

これ、ひょっとしてCreateCliqueTreeの中でeliminateされるんじゃね？
とりあえずそう信じて書いて動かしてみよう。

＞違った

とりあえず前回のDesicionMakingの課題のPrintFactorとかを持ってきてプリントしたりしてみよう。
90：karino2：2017/01/09(月) 22:14:29: ネットワークの図を眺めながら思った事。
factorList(2)はP(2|1, 3)だよな。
一方でクラスターグラフの2, 8はscopeが2, 8のファクターを入れるのが正しいよな。
2, 8のスコープのファクターといえばP(8 | 2)だよなぁ。

そう考えるとこのcliqueに入るのってP(8 | 2)だけな気がするんだが。
と思って眺めていると、RESULTとはvarの順番が違うだけな事に気づく。

ほぉ。なるほど。

つまり、ファクターのアサインは自分でやらなきゃいけない、という事じゃないか。
どうアサインしたら良いのだろう？とりあえず先頭から何も考えずにアサインしてみよう。
これだと全部のクラスタにアサインされる、という風にはならん気もするが、そういう風にクラスタが作られてると思ってみよう（自分には理解出来ていないが）。

という事でpdfのアルファがどうこう、と言っている意味をようやく理解する。
FactorProductを参考にマップとか作って解く。
無事とけた。なんかここまでで3時間くらいかかってるな…
91：karino2：2017/01/09(月) 22:17:56: GetNextCliques.m を実装する。
まずサンプルを見て途方に暮れる。
どうもINPUT1を第一引数に、INPUT2を第二引数にする、という事か？

で、適当にINPUT2を見るがだいたい空。ふむ。
良く分からないのでGetNextCliques.mのコメントを読む。

読んだ。どういうメッセージが来うるか、という事をどうしるんだろう？edgesか。
ふむ、

という事はあるクリークiとつながりのあるクリークはedges(i, :)で分かるのかな。
そんな感じで実装してみよう。

何も考えずに実装したら7, 4となったが、サンプルは9, 2となっている。
そうか。送った物じゃなくて受け取った物を考えないといけないのか。

なんかメッセージを一つ受け取ってない場合と全部受け取ってる場合の二つのif文が出来てしまったが、とりあえず動いている模様。
92：karino2：2017/01/09(月) 23:14:48: ComputeExactMarginalsBPまで終わった！
大変過ぎる…
93：karino2：2017/01/09(月) 23:38:55: Max-Sumやろうとしたが、これMAP推定の動画見終わってからの方が良さそうだな。
という事でここまででひとまずおしまい！
94：karino2：2017/01/10(火) 11:03:03: 見終わって実装してみた。こちらは簡単ね。
あっさり最後まで終わる。やったね！
95：karino2：2017/01/12(木) 07:19:02: 明け方に目が覚めてしまったので、暇つぶしにサンプリングの動画を見ていた。
Gibbsサンプラーは理解。
で、Metropolis-Heisting法の途中までだが、詳細釣り合い条件が凄い良くわかった。
このコースのMCMC関連の講義動画の出来は感動的だな。
96：karino2：2017/01/12(木) 10:36:48: MHアルゴリズムの問題も見終わった。
ちゃんと理解出来たと思うが、これはプログラム課題やってみたくなるね！

あとでやろう。
97：karino2：2017/01/20(金) 17:30:54: MCMCのPAやるです。
開幕のBlockLogDistributionの入力が分からない。
exampleINPUTのt1a1とt1a2じゃないの？と思うが、全然関係ないように見える。
そしてTestToyはComputeExactMargjnalsBPが無くて途中から実行出来ない。前回のPAでやったが何を持ってきたらいいか分からんのでとりあえずおいとく。

TestToy.mでtoy_networkとtoy_factorとtoy_evidenceは作れるので、これとsubmit.mをみくらべると、[1]をVに入れると呼べそうだな。
期待値は分からんが。
98：karino2：2017/01/20(金) 17:35:19: exampleINPUTを眺めていたら、t6a2がグラフっぽいな。
なんかt6a1からt6a4までがちょうど引数になってるように見える。
exampleOUTPUTのt6が結果か？
なんか答えが2つあるが。サンプリングなら変数一つにつき結果一つじゃないの？

とりあえず何をやる関数か考えてみよう。
99：karino2：2017/01/20(金) 17:49:02: 基本的には、全ファクターにevidenceと一致するセルを集めてグラフ構造に従って掛けて、Vの分布テーブルを作って返す、という事じゃないか。
ノーマライズされてないのだから、Vが含まれていないファクターは関係ないか。
まずマルコフブランケットをGからどう取ってくるか、を解読するのかな。
100：karino2：2017/01/20(金) 18:21:58: var2factors がそれっぽい。
Vはベクトルだが、結果は全部の変数に同じ数が入ってるのを仮定してよさそう。
Vの入ってる全ファクター取ってきて、V以外の要素をevidence でルックアップすると、Vのファクターとなる。
ここから全部が同じ値の物だけ取っていってlogすれば良い、、、かな？
101：karino2：2017/01/20(金) 20:53:37: メモいろいろ。

>> a = repmat((1:3)', 1, 2)
a =

1 1
2 2
3 3

>> [a, ones(3, 1)]
ans =

1 1 1
2 2 1
3 3 1

>> b = [2, 3, 4]
b =

2 3 4

>> b .* ones(3, 1)
ans =

2 3 4
2 3 4
2 3 4

>> a = zeros(5, 3)
a =

0 0 0
0 0 0
0 0 0
0 0 0
0 0 0

>> a(:, 2) = ones(5, 1)
a =

0 1 0
0 1 0
0 1 0
0 1 0
0 1 0

>> a(:, [2, 3]) = ones(5, 2)
a =

0 1 1
0 1 1
0 1 1
0 1 1
0 1 1

>> BlockLogDistribution(ex.t6a1, ex.t6a2, ex.t6a3, ex.t6a4)
ans =

0.00000 0.40547

>> exampleOUTPUT.t6
ans =

0.00000 0.40547

お、一致した。これで良さそうね。
102：karino2：2017/02/07(火) 13:11:54: GibbsTransは randi('seed',1) した後やると結果が一致。

>> load("exampleIOPA5");
>> ex = exampleINPUT;
>> eo = exampleOUTPUT ;

>> [res1, res2] = MCMCInference(ex.t8a1{1}, ex.t8a2{1}, ex.t8a3{1}, ex.t8a4{1}, ex.t8a5{1}, ex.t8a6{1}, ex.t8a7{1}, ex.t8a8{1});

よし、次はMHUniformTrans、と思ったがこれが全然分からない。
定常分布がなんなのか分からない。
マルコフネットワークの定常分布とはなんぞや？

マルコフネットワークの定義に戻ると、全ポテンシャルをかけた物が同時確率となる訳だ。
これから定常分布を求めるとはどういう事だろう？
いや、これこそが定常分布だよな。
それに対して提案分布があって、提案分布でサンプリングしたランダム変数の組みが同時確率としてどれくらい成立しうるか、という話じゃないか？

お、わかった気がする。

res4 = MHUniformTrans(ex.t9a1{4}, ex.t9a2{4}, ex.t9a3{4});
eo.t9{4}

うーん、毎回結果が似すぎていて良く分からないな。
103：karino2：2017/02/07(火) 13:18:33: レシート入力を終えて久しぶりに再開。

>> load("exampleIOPA5");
>> ex = exampleINPUT;
>> eo = exampleOUTPUT ;

res4 = MHUniformTrans(ex.t9a1{4}, ex.t9a2{4}, ex.t9a3{4});
res4
eo.t9{4}

なんか合ってる気がする。とりあえずMHUniformTransを読み直してみよう。

思い出してきた。。Qはユニフォームな分布なので、x->x'ももx'->xも同じ確率。
で、pi_xはfの積でassignmentはAとA_prop。
うむ、この式であってそうだな。

MHSWTrans.mの実装に入る。
pdfの図をみるとアルゴリズムの概要は分かる。
ただ、具体的に実装しようとすると難しい。

実行してみるとscomponentsが無い。そこで何が必要か逆算で考えてみる。

>> SelEdgeMat
SelEdgeMat =

Compressed Column Sparse (rows = 16, cols = 16, nnz = 14 [5.5%])

(4, 3) -> 1
(7, 3) -> 1
(3, 4) -> 1
(10, 6) -> 1
(3, 7) -> 1
(10, 9) -> 1
(6, 10) -> 1
(9, 10) -> 1
(15, 11) -> 1
(16, 12) -> 1
(11, 15) -> 1
(16, 15) -> 1
(12, 16) -> 1
(15, 16) -> 1

ここから
[var2comp, cc_sizes] = scomponents(SelEdgeMat);
が欲しい。

ccはあとのコードを読むとconnected componentの略みたい。
var2compoは

selected_vars = find(var2comp == selected_cc);

ででconnected componentの変数一覧が取れるようなもの。
cc_sizesは各connected componentのサイズが入っているのだろう。
var2compはvarと同じ長さのベクトルで値は属しているccのインデックスだろう。

selected_edgeはpdfで言う所のactivatedという奴だな。

さて、SelEdgeMatからCCはどうやったら出るだろう？

1から順番に見ていき、繋がってるエッジを辿ってそいつをCC1にアサイン、繋がってたノードはスキップして2に進む、という感じか。

scomponentsの実装は終わり、、MHSWTransも実装してみたが、なんかこれ、サンプルが毎回入力と出力一緒だなぁ。

>> MHSWTrans(ex.t10a1{1}, ex.t10a2{1}, ex.t10a3{1}, ex.t10a4{1})

pairwiseのfactorを考慮に入れた遷移がちょっと面倒なので、現実逃避にDiscussion Forumなどに答えたりscomponentsの自分実装を貼ったりする。
104：karino2：2017/02/07(火) 13:19:34: scomponentsはgaimc下にあるぜ、との事。
addpath('gaimc');
で良さそう。せっかく書いたのに(´・ω・｀)

さて、ちょっと確定申告やってたら、ここまでの話を忘れた。どうしよう？
まずGってなんだっけ？

edgesにi, jにエッジがあれば1、そうでなければ0が入っているのだろう。

var2factorsってなんだっけ？ >>100 に何か言及があるが。

>> ex.t10a2{1}.var2factors{1}
ans =

1 17 20

変数1が含まれているファクターは1, 17, 20という事か？

ex.t10a3{1}(1)
ex.t10a3{1}(17)
ex.t10a3{1}(20)
などを試してそうっぽい、と納得。

q_listってなんだろう？
一列目と二列目で、エッジのある頂点の組みを表している。
三番目の値はランダムに切る時に使われているだけなので、意味は無いんじゃないか（全部0.5だし）。

MCMCInferenceを見て分かった。q_listは自分で作るものだ。q_ijだね。
なるほど。

-----

BlockLogDistributionの実装が不完全で直したり。

>> randi('seed',1);
>> MHSWTrans(ex.t11a1{2}, ex.t11a2{2}, ex.t11a3{2}, ex.t11a4{2})
ans =

2 2 2 2 1 2 2 1 2 2 1 1 1 2 1 1

>> eo.t11{2}
ans =

2 2 2 2 1 2 2 1 2 2 1 1 1 2 1 1

治ったかな？
105：karino2：2017/02/07(火) 13:20:37: MCMCInferenceのq_ijに進む。まずはEdgeToFactorCorrespondence(G.names, F)が何を返すか調べるか。
>> E2F = EdgeToFactorCorrespondence(ex.t12a1{1}.names, ex.t12a2{1})
E2F =
{
[1,1] = [](0x0)
[2,1] = 17
[3,1] = [](0x0)
[4,1] = [](0x0)
...

つまりエッジに対して、そのエッジを含むfactorを返すみたいね。

なんとなく実装して試すと、exampleINPUTの1のMHUniformの結果がそもそもあってない。あれ？これって前実装した処だよなぁ、とC-rしてみたらex.t8a1とかで呼び出している。
もう一回見てみるか。

[res1, res2] = MCMCInference(ex.t8a1{1}, ex.t8a2{1}, ex.t8a3{1}, ex.t8a4{1}, ex.t8a5{1}, ex.t8a6{1}, ex.t8a7{1}, ex.t8a8{1});

＞うーん、乱数のseed設定したらex.t8aの結果とeo.t8o1は厳密に一致しているなぁ。

ex.t12a*{1}とex.t12a*{2}がeo.t12o1{1}とeo.t12o1{2}に一致してないが、やっぱり{1}が一致しないという事が良く分からない。
UniformTransのコードを見直したが、正直自分のところにバグが入る余地は無い気がするので、検証の仕方が間違いなんじゃないか、という事で先に進んでみる。

----

TestToy.mの下の方のMCMCInferenceのブロックを実行したら帰ってこない。

凄い経ったら帰ってきたが、ExactMが無い、とか言って途中で終わってしまった。ComputeExactMarginalsBPの結果か。
でもこれ、時間かかり過ぎて検証出来ないなぁ。
どうしたもんか。
106：karino2：2017/02/07(火) 21:02:17: やはりこれは自分のマシンでは動作確認は諦めよう。
理屈は理解したので、この実装が例え間違っていても、実際に実装しなきゃいけない時には実装出来ると思う。
107：karino2：2017/02/10(金) 14:28:38: dual decompositionとtractable MAP inferenceの動画を見てなかったので観ている。
これが終わればコース2は終わりでいいかな。