予備実験のために, 単純なディリクレ分布(DMのようにmixtureでない もの)の推定の計算をする必要があった。
1年位前に書いたコードがあるのだが, 先日の疎行列のフォーマットを 使ってコードを書き直してみた。

DMでは LOO (Leave-one-out) 近似を使っているが, その際気になったので, Ψ関数を使った exact な計算と, MacKayの近似法(MacKay 1994. 去年の SVM2004で紹介した話)を比べてみた。
小さなデータの代表として Cranfield コーパスを使った結果, 下のような感じになった。

	Unigram	Bigram	Trigram
Exact	- (49860)	64.20 [31] (57.2312)	- (8298.3)
MacKay	385.99 [15] (240.1676)	236.80 [9] (45.0992)	416.09 [16] (195.6051)
LOO	24.30 [40] (293.6186)	14.17 [24] (57.5734)	29.55 [34] (24.3058)

数値は順に, 計算時間(秒), iterationの回数, α = Σ_kα_k. - は収束しなかったことを表す。
UnigramはDMで混合数 M=1 の場合, Bigramは (MacKay 1994) の場合と同じ。

これを見てわかるように, bigram以外では, Exact な iteration は収束しない。 山本先生のところ で, DMの計算にΨ関数を使った Exact な計算では収束しなかった, とおっしゃって いた記憶があるが, DMでなくて単純なディリクレ分布の場合でも, Exactな fixed-point iteration は unigram については収束しないようだ。
(()の中は100回の最大の繰り返し終了後の値だが, 見ると, αがすごい値になってしまっていることがわかる。)

最適化の振る舞いを見ていると, 一回最適解に行きかけるが, 途中でそれを超えて 発散してしまっているように見える。理論的には単峰で凸なはずなので 謎。ただ, 他のデータでもそうだったので, 山本先生のグループのおっしゃっていた ことは正しいのだと思う。
という, Minkaの論文を見ていただけではわからない実際的な話がわかったので, まとめてみました。

プログラムは間違っていないと思うのだが, 参考のために置いておきます。
(あくまで参考なので, これだけでは動かないので注意。)
dirpolya.m dirpolya_mackay.m dirpolya_loo.m

fastfit

MacKay本 の612ページをたまたま見ていたら,

e⁷ = 1096

ということを知った。へー。 つまり, log(1100) 〜 7 くらい, ということ。(底は e)

ちなみに, 正確には

• log(1100) = 7.0030
• log(1000) = 6.9077

くらいのようだ。つまり, log(1000) 〜 6.91 くらい。 (exp(6.91) = 1002.2)

何が嬉しいかというと, log(1100) 〜 7 を参考にして, もう少し小さい数の log(1000) 〜 6.91 を覚えておけば,
log(10^3) = 3log(10) 〜 6.91 ∴ log(10) 〜 6.91/3 = 2.3 (2.303) なので,
log(10000)は何? と言われたら, log(10000) = log(10^4) = 4log(10) 〜 9.2 くらい, と簡単にわかる ということ。 つまり, 10^x くらいの数の自然対数は, 2.3x くらいになるとわかる。
逆に, log(x) = y がわかっているとき, この x は 10^(y/2) くらい, ということ。 計算機に入っている log はたいてい自然対数なので, この関係は 役に立つと思う。たとえば, x の確率が log(p(x)) = -1000 だったとすると, 実際は p(x) ~ 10^(-500) くらいかなとわかる。

log(x) = log₂(x)/log₂(10)なので, log₂(10) = 3.321を覚えておけば, 2^10 = 1024 を使うという方法もあるが, もっと直接的な 方法ということで。
「ご冗談でしょう、ファインマンさん」とかを読んでいると, 物理の人にはこういうのは常識なのかも知れないと思うけれども。

明日7/22の DMLA では, Unsupervised Boosting の話をします。
具体的には, 以下の論文を読みます。
Max Welling, Richard S. Zemel, and Geoffrey E. Hinton. "Self Supervised Boosting". NIPS 2002. [LINK]

興味のある方はご参加ください。NAIST A707, 17:00- です。

ぐはっ。難しい..。
でも, だいぶわかってきた。
以下, 公開デバッグというか, 自分用メモ。

• 確率変数 X とは, 通常は Borel 集合からRへの関数。 (Borel集合全体の measure は 1)
• 確率変数 X に対して, 確率分布 P_X(B) とは, 値域 B をいろいろ変えたときの, 定義域 X^-1(B)のルベーグ測度を返す関数 (→[0,1]). だから, P_X(B)も1つの測度。
• Xがconstant: X≡a であれば, P_X(B)は Dirac測度 δ_a(B)になる。
• ルベーグ測度は, ルベーグ外測度を使って定義できる。
• 関数 f が Lebesgue 可測 (measurable,ルベーグ積分可能) とは, fの値域に属する任意の interval I ⊆ R に対して, その逆像 f^-1(I) がルベーグ可測であること。
• 集合がルベーグ可測とは, 要するに, そのクラスがσ-fieldになっている, ということ。(普通の interval はみんなルベーグ可測)

特に X が二値で, X = a (probability p), X = b (probability 1-p) であるとき, 確率分布 P_X(B) は p δ_a(B) + (1-p) δ_b(B)に当然なる。
一般に, 確率変数(=全体の測度が1であるBorel集合からRへの関数) X の 値が離散的に, a_i を確率 p_i でとるとき (i = 1,2,..),

P_X(B) = Σ_i=1^∞ p_i δ_{a_i}(B)

と自然にかける。

ここで, 引数Bとなる, Xによって生成されるσ-fieldの要素は, Xの解像度によって制限があることに注意。Xが {a,b} への関数であれば, Xによって生成される σ-field は {0,Ω,X^-1({a}),X^-1({b})}になる。
ってところなのかな。