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

• gamglm, Bayesian Gamma generalized linear model.

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

• 特徴ベクトル f(n)=(f⁽ⁿ⁾₁, f⁽ⁿ⁾₂, ..., f⁽ⁿ⁾_M) に対して, 重みベクトル w_aとw_bを用いて
  • a(n) = exp(w_a^Tf(n)), b(n) = exp(w_b^Tf(n)),
• 観測値 y(n) ~ Ga(a(n),b(n)).

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

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

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のようなデータで正の実数を回帰するプログラムとして, 使われることがあれば よいなと思っています。

---

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

相変わらず, Γ関数の海を泳いでいるわけですが..。
いま, 語彙が全部で V 個あったとして, ある文書において, それぞれの単語が x = (n_1, n_2, ..., n_V) 回現れたとする。 単語のユニグラム確率を p_i とすれば, このデータ x の確率は p_i が n_i 回現れているので, 多項分布に従って

p(x|p) = (n₀!/Π_in_i!)Π_i p_i^n_i　　　　　(1)

になる。ただし, n₀ = Σ_i n_i.
これを最大にする確率 p_i は最尤推定で,

p_i = n_i/n₀　　　　　(2)

として求まる。 これは基本的に Naive Bayes で使われている方法で, 高村君のNBのチュートリアル などを参照。
ただ, n_i はほとんどが0で非零要素はせいぜい100個くらいなのに, 10000次元ぐらいある p_i が正確に求まるはずがない。 つまり, 確率 p_i 自体に不確実性があるということなので, p_i を生成する 確率分布 p(p|α) をもってきて,

p(x|α) = ∫p(x|p)p(p|α)dp

として p を積分消去してしまうことができる。この積分は解析的に計算できて,

p(x|α) = Γ(α₀)/Γ(α₀+n₀)Π_i Γ(α_i+n_i)/Γ(α_i)　　　　　(4)

になる。この分布は Polya分布 (Dirichlet/Multinomial)と呼ばれている。 ただし, α₀ = Σ_iα_i.
ここまでが前フリ。

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

いま, (1)式は

p(x|p) ∝ Π_ip_i^n_i

を意味するので, 両辺の log をとると

logp(x|p) = Σ_i n_i log p_i.

よって, ∂logp(x|p)/∂p_i = n_i/p_i.つまり, n_i = (∂log(x|p)/∂p_i)・p_iということになる。
これを Polya 分布の場合に同じ計算をしてみると,

(∂logp(x|α)/∂m_i)・m_i = α_i[Ψ(α_i+n_i) - Ψ(α_i)] ≡ n~(i).

になる。ここで, m_i = α_i/α₀.
これは n_i と同じような意味を持っていると考えられるので, n~(i) とおくと, これがどういう意味を持っているかが問題になります。 以下の話は, Minkaの論文には書いていない話。

さて, digamma関数 Ψ(x) = d/dx logΓ(x) は, Ψ(x+1) = 1/x + Ψ(x) という性質を持っているので (Γ(x+1) = xΓ(x) なので, 両辺を対数微分すれば簡単に出てきます。 きのうの夜はこれに気付かずに, ワイエルシュトラスの公式を使って一生懸命計算して しまった。; ), 実は Ψ(x+n) - Ψ(x) = Σ_i=1ⁿ 1/(x+n-i)です。つまり, n~(i) は (n_i = n, α_i = αと書くと)

n~(i) = α Σ_i=1ⁿ 1/(α+n-i)
= Σ_i=1ⁿ α/(α+n-i)
= Σ_i=1ⁿ 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 が小さい→ x_i の事前観測値が小さければ, x_i が何回観測されても, 1回だけ観測されたのと同じになる, ということ。 これが Minka の論文の(79)式の意味する所だと思う(こういう説明は論文には 全然書いていないが)。

LOO approximation

MacKay's interpretation