mots quotidiens.
Daichi Mochihashi (持橋大地) daichi <at> ism.ac.jp by hns, version 2.10-pl1.

先月 2024年04月 来月
1 2 3 4 5 6
7 8 9 10 11 12 13
14 15 16 17 18 19 20
21 22 23 24 25 26 27
28 29 30

2014年10月28日(火) [n年日記]

#1 ガンマ一般化線形モデル

共同研究の関係で必要になったので, ガンマ分布の一般化線形モデルのC++ソフトウェア を公開しました。

ガンマ一般化線形モデルは, 以下のようなモデルです。

一般化線形モデル自体の説明については, 大慈さんの講義スライド などがわかりやすいと思います。

普通, 言語処理などでは次のようなデータファイルを使って

1  w0|の p0|助詞 w-1:家 p-1|名詞
0  w0|人 p0|名詞 w-1:の p-1|助詞
先頭の二値の観測値をSVMやロジスティック回帰で予測します。しかし, この観測値が 0/1ではなく, 0以上の実数である場合も普通にありえます。
この場合, 観測値がガンマ分布 Ga(a,b) に従うとするのが自然なモデル化の一つで, a > 0, b > 0であるため, ロジスティック回帰のように a = exp(<w,f>), b = exp(<v,f>) のように指数の肩に乗せてやることができます。 *1 上の gamglm はこのためのソフトウェアで,
y  feature_1 feature_2 feature_3 .. feature_n
のようなデータファイルに対して動きます。yは任意の正の実数です。

ガンマ一般化線形モデルはRの glm() で推定できるような気がしてしまいますが, パラメータ w に対して凸ではないため, Webで例にあるような説明変数の数が たかだか1個か2個の場合ではなく, 何百, 何千とある場合はまともに動きません (のはず)。実際, 微分を計算して L-BFGS で最適化するプログラムを最初書いたところ, すぐ局所解にはまってしまい, 変数の数が増えると結果がほとんど0になるなど, 使いものにならない結果となることがわかりました。

このため, 上の gamglm では, 地道なランダムウォークMCMCでパラメータを推定して います。ただし, 特徴の数が数千を超えると尤度の計算が非常に遅くなるため, 特徴が多くの場合疎であることを利用して, 現在更新する特徴が現れているデータ だけを利用してMetropolis-Hastingsの尤度計算をする, ということをしています。
やっていることは簡単ですが, 効率的に実装するには結構時間もかかり, 一般的な問題でもあるので, 公開することにしました。
SVMのようなデータで正の実数を回帰するプログラムとして, 使われることがあれば よいなと思っています。


*1: 一般にR等でガンマ線形モデルのリンク関数は逆数 1/x にすることが多いようですが, これは正値であることを何ら保証しないので, 変量の数が1や2ではない場合は使いものにはならない気がします。

2005年10月28日(金) [n年日記]

#1 Polya分布再考

下に色々書き足したので, ここまで来ると, 別に文書にした方がいいような 気がしてきました。
と言っても, そんなに完全に新しいことを言っているわけ ではないと思うので, 少し躊躇があります。

相変わらず, Γ関数の海を泳いでいるわけですが..。
いま, 語彙が全部で V 個あったとして, ある文書において, それぞれの単語が x = (n_1, n_2, ..., n_V) 回現れたとする。 単語のユニグラム確率を p_i とすれば, このデータ x の確率は p_i が n_i 回現れているので, 多項分布に従って
p(x|p) = (n0!/Πini!)Πi pi^ni     (1)
になる。ただし, n0 = Σi ni.
これを最大にする確率 p_i は最尤推定で,
pi = ni/n0     (2)
として求まる。 これは基本的に Naive Bayes で使われている方法で, 高村君のNBのチュートリアル などを参照。
ただ, ni はほとんどが0で非零要素はせいぜい100個くらいなのに, 10000次元ぐらいある p_i が正確に求まるはずがない。 つまり, 確率 p_i 自体に不確実性があるということなので, p_i を生成する 確率分布 p(p|α) をもってきて,
p(x|α) = ∫p(x|p)p(p|α)dp
として p を積分消去してしまうことができる。この積分は解析的に計算できて,
p(x|α) = Γ(α0)/Γ(α0+n0i Γ(αi+ni)/Γ(αi)     (4)
になる。この分布は Polya分布 (Dirichlet/Multinomial)と呼ばれている。 ただし, α0 = Σiαi.
ここまでが前フリ。

Γ関数を含むこの分布は僕は毎日のように使っているわけだけれども, 一体この分布はどういう意味を持っているのだろうか, というのが気になった。
実はこの答えは, Minka の "Estimating a Dirichlet distribution" の4章に書いてあるのだが, とりあえず必要がなかったので, これまで 読み飛ばしていた。(この論文は普通に見るとかなり難しいと思うので, 隅から隅まで読んで計算をフォローしている人はあまりいないと思うのですが..。) ここに書いてある議論は, 以下のようなこと。

いま, (1)式は

p(x|p) ∝ Πipi^ni
を意味するので, 両辺の log をとると
logp(x|p) = Σi ni log pi.
よって, ∂logp(x|p)/∂pi = ni/pi.つまり, ni = (∂log(x|p)/∂pi)・piということになる。
これを Polya 分布の場合に同じ計算をしてみると,
(∂logp(x|α)/∂mi)・mi = αi[Ψ(αi+ni) - Ψ(αi)] ≡ n~(i).
になる。ここで, mi = αi0.
これは ni と同じような意味を持っていると考えられるので, n~(i) とおくと, これがどういう意味を持っているかが問題になります。 以下の話は, Minkaの論文には書いていない話。

さて, digamma関数 Ψ(x) = d/dx logΓ(x) は, Ψ(x+1) = 1/x + Ψ(x) という性質を持っているので (Γ(x+1) = xΓ(x) なので, 両辺を対数微分すれば簡単に出てきます。 きのうの夜はこれに気付かずに, ワイエルシュトラスの公式を使って一生懸命計算して しまった。; ), 実は Ψ(x+n) - Ψ(x) = Σi=1n 1/(x+n-i)です。つまり, n~(i) は (n_i = n, α_i = αと書くと)
n~(i) = α Σi=1n 1/(α+n-i)
= Σi=1n α/(α+n-i)
= Σi=1n 1/(1+(n-i)/α).
ここで, Σの中の各項が1に等しければ, 当然 n~(i) = n_i です。 しかし実際は, 各項は
1/(1+(n-i)/α)
で, 1より小さくなっています。つまり, この和は
1 + 1/(1+1/α) + 1/(1+2/α) + ... + 1/(1+(n-1)/α)
になっていて, 頻度 i が大きくなるほど, その実際のカウントは1増えるのでは なく, 1/(1+(n-i)/α) = 1 - (n-i)/(α+n-i) だけ "dumping" されます。 Minka の論文の Figure 1 には, この和がαが小さい時には log のような ダンピング効果を持つと書いてありますが ∫1/xdx = log(x) なことを考えると (この場合は離散ですが), わりと納得できると思う。
つまり, ある単語が20回出たという事象は1回出た事象の20倍の意味がある (普通の多項モデル)ではなく, log(20)倍くらいだろう (対数の底は任意として), という常識的な直観をモデル化していることになっている。 もちろん, この dumping は αi, つまり事前の仮想的な観測数に依存していて, αiが大きい, つまり観測が多いと仮定されていればほとんどダンピングはなく 線形に伸びるが, αiが小さく観測があまり信頼できないと予想されていれば, たとえその後何回も観測された としても, その効果は log 程度に緩められる, ということになる。

実際, ここで α が 0.001 のように非常に小さい時(自然言語の場合は よくある), 上の級数の2項めから後ろは分母が非常に大きくなるのでほとんど0に 近付いて, 結局「n が何であっても」この effective count はほぼ1になる, ということ。(!)
つまり, もし αi が小さい→ xi の事前観測値が小さければ, xi が何回観測されても, 1回だけ観測されたのと同じになる, ということ。 これが Minka の論文の(79)式の意味する所だと思う(こういう説明は論文には 全然書いていないが)。

・ LOO approximation

文書データ X = (x1, x2, .., xN)が与えられた時, 上の(4)式のαを最適化するためには,
αk' = αk ・(ΣiΨ(αk+nik)-Ψ(αk)) /(ΣiΨ(α0+ni)-Ψ(α0))
という fixed-point iteration を計算します。 (この式を導くには, 色々自明でない bound を使って計算する必要があります。)
ただ, DM の論文にあるように, この計算はΨ関数を使っていて重いので, LOO近似を使うと, また自明でない bound を使って結構計算して,
αk' = αk ・(Σi nik/(αk+nik-1)) /(Σi ni/(α0+ni-1))
という fixed-point iteration を解けばいいことがわかります。
この式の意味は僕はずっとよくわからなかったのですが, 上で書いたように, Ψ(x+n) - Ψ(x) = Σi=1n 1/(x+n-i)なので, 実はこの近似は
Ψ(α+n) - Ψ(α) = Σi=1n 1/(α+n-i) ~ n/(α+n-1)
という近似をやっているのだ, ということに気がつきました。
つまり, 1/(α+n-1) + 1/(α+n-2) + .. + 1/α という和を, 各項をみな一番小さい値 の 1/(α+n-1) で近似して, n * 1/(α+n-1) として計算しているということ。 これだと誤差が出るはずですが, 分子でも分母でも同じ近似をしているので相殺されて いるのかもしれません。この近似が, LOO近似は
Γ(x+n)/Γ(x) = x(n) = x(x+1)..(x+n-1) ~ (x+n-1)^n
という近似をしているのだ, という Minka の(71)式の意味なのだと思います。 (この式には説明がないので, 前に読んだ時は(71)式の理由が全然わからなかった。)

・ MacKay's interpretation

実は MacKay も, "Models for Dice Factories and Amino Acid Probabilities" (PDF)の最後で少し違った解釈を書いています。
ただ, MacKayの近似は α_i > 0.1 のような大きい場合を対象にして計算して いますが, 言語の場合は α_i はずっと小さいのであまり当てはまらないような 気はします。小さい場合の近似式を使って計算してみたところ, あまり 面白くない式になってしまいました。(大きい場合には, KLダイバージェンス*αで 正則化がかかっていると解釈できるような式が得られる。)


2 days displayed.
タイトル一覧
カテゴリ分類
 なかのひと
Powered by hns-2.10-pl1, HyperNikkiSystem Project