とある童話の統計科学(3)言語モデルとn-gram

Toaru_Logo

そういえば、この辺の話は基本的に、C. Manning, “Foundation of Statistical Natural Language Processing”に詳しい。ぜひご参考に。

言語モデル

言語研究は統計に限らないので、定義はまちまちだろうけども、ここでの言語モデルとは、語の有限列である文書$W = \{w_1, w_2, \cdots, w_N \}$に対して、確率$p(W)$を与える関数を指す。ただし、語はボキャブラリ数をVとして、自然数にコーディングされているものとする、すなわち、$w_i \in \{1, \cdots, V \}$である。

例えば、次のような2つの英文があったとしよう。

  • an apple is red
  • in are take me

1つ目は、単に「りんごは赤い」と言っているだけである。一方で、2つ目の文章は何を言いたいのかさっぱりわからない。「英語」という言語モデルにとって、p(an apple is red)はある程度の確率を持つが、p(in are take me)の確率は極めて小さいとも言える。

統計的言語処理の1つの目的は、膨大に蓄積された言語リソースから、この言語モデルを統計的に推定することと言える。今回のテーマであるn-gramとは、その言語モデルの1つである。

n-gram

n-gramが推定するモデルは、 p(w_i | w_{i-1}, \cdots, w_{i-n+1}) である。すなわち、直近のn-1語によって、次に出現する語の確率を与える。例えば、次の2つの文章を考えてみる。

  • this is ***
  • he eats ***

***に何が入るか、予測できるだろうか?1文目は、これは〜です、と言っているため、候補は非常に多い。一方、2文目は彼は〜を食べます、なので〜に入るものは食べ物と予測できる、かつ「彼」という言葉からは男性に一般的に好かれるような食べ物が入る可能性が高い。この我々の中に内在する「直前の語からなんとなく予測できる」が、n-gramで抽出したい言語の特性である。

この問題は形式的には、「直前2語を所与とした場合の、次に出現する語の条件つき確率を求めよ」、ということである。

  • p(*** | this is)
  • p(*** | he eats)

この例は、直近の2語から次に出る1語を予測する問題なので、3-gram(trigram)と呼ぶ。

n-gramに対して、文書Wの確率は、次のように与えられる。

(1)    \begin{align*} p(W) = \prod_i p(w_i | w_{i-1}, \cdots, w_{i-n+1}) \end{align*}

また、特殊な形として、1-gram(unigram)がある。これは、直前の語による影響を受けないモデルである。この場合の文書の確率は以下である。

(2)    \begin{align*} p(W) = \prod_i p(w_i) \end{align*}

そのうち紹介するLDAなどのトピックモデルは、この1-gramの発展系である。

n-gramの推定

簡単のために、2-gramを例とする。ベイズの定理より、次が成り立つ。

(3)    \begin{align*} p(w_n | w_{n-1}) = \frac{p(w_n, w_{n-1})}{\sum_{w_n} p(w_n, w_{n-1})} \end{align*}

ここで、$p(w_n, w_{n-1})$として、パラメータ$\{ \theta_{i,j} \}$、観測数1回の多項分布を考える。すなわち、

(4)    \begin{align*} p(w_n=i, w_{n-1}=j) = \theta_{i,j}, \sum_{i,j} \theta_{i,j} = 1, \theta{i,j} > 0  \end{align*}

文書Wに対して、パラメータ$\{ \theta_{i,j} \}$の最尤推定量は、

(5)    \begin{align*} \theta_{i,j} = n_{i,j} / N \end{align*}

$n_{i,j}$は、語列{i,j}の出現回数、Nはドキュメント内の総語数である。ここで、Eq.(3)の2-gramモデルによる条件つき確率は、次のように計算できる。

(6)    \begin{align*} p(w_n=i | w_{n-1}=j) = \frac{\theta_{i,j}}{\sum_{\tilde{i}} \theta_{\tilde{i},j}} = \frac{n_{i,j}}{\sum_{\tilde{i}} n_{\tilde{i},j}} \end{align*}

適用例

以下が、2-gramのパラメータ推定をしたのちに、各文書についてp(*|she), p(*|red), p(*|look)の確率が高い上位5語をそれぞれ抽出した結果である。

ngram_ex

主語となるsheの後には動詞が続いていること、形容詞redの次には赤いもの、動詞lookの次にはatなどの前置詞が続き、我々の英語知識に沿った結果が得られているのではないだろうか。

今回の2-gramのパラメータ総数はV*V=121,396,324個である。しかし、グリム童話の場合、そのうちの実に99.97%のパラメータが値は0のまま、つまりその組み合わせが本文中に観測できなかった。このスパースさが、n-gramの難しい所である。概念的には同じようなものを指していても、最尤推定されたパラメータによる文書の確率推定は度々0になり(Eq.(1)のうち1つでも0になれば、総積によって全体が0となる)、モデルとして役に立たない。文書の量を増やしていくというのは1つの手段ではあるが、3-gram, 4-gramを考えていったとき、パラメータの数は指数的に増大していくため、本質的な解決にはならない。

次回は、その緩和としてスムージングと呼ばれる手法の紹介をする。

無限混合正規分布 cont.

前回の補足。

(たぶん)実装可能なメモ書きがこちら

psuedoコード代わりのRコード。
データを入れれば、一応動くけれど、ナイーブな実装なので超遅い(Rでも)。

# 1. prepare dataset
X = ...

# 2. set up hyperparameters
theta0 = list(NULL)
theta0$m = c(0,0)
theta0$xi = 1e-2
theta0$eta = 2*D+2+1e-1
theta0$B = 1e-1*diag(D)*1

# 3. Gibbs sampling
DPM = function(X, theta0){

	# 3.1 functions

	# logarithm of multivariate gamma function
	ln_mgamma = function(x, d){
		0.25 * d * (d - 1) * log(pi) + sum(lgamma(x + 0.5 * (1 - 1:d)))
	}

	# calculate posterior of NIW
	posterior = function(X, theta0){

		if(class(X) == "numeric"){
			X = t(X)
		}

		N = nrow(X)
		D = ncol(X)

		theta = list(NULL)
		theta$xi = theta0$xi + N
		theta$m = (theta0$xi * theta0$m + apply(X, 2, sum)) / theta$xi
		theta$eta = theta0$eta + N
		theta$B = theta0$B + theta0$xi * (theta0$m %*% t(theta0$m)) + t(X) %*% X - theta$xi * (theta$m %*% t(theta$m))

		return(theta)

	}

	# logarithm of the gaussian marginal likelihood for NIW parameter distribution
	ln_marginal_likelihood = function(X, theta0){

		if(class(X) == "numeric"){
			X = t(X)
		}

		N = nrow(X)
		D = ncol(X)

		theta = posterior(X, theta0)

		a = 0.5 * (theta0$eta - D - 1)
		ap = 0.5 * (theta$eta - D - 1)
	
		retval = -0.5*D*N*log(pi) + 0.5*D*log(theta0$xi/theta$xi) + ln_mgamma(ap, D) - ln_mgamma(a, D) + 
				a * determinant(theta0$B)$modulus - ap * determinant(theta$B)$modulus

		return(retval)

	}

	# 3-2. initialize MCMC state & parameters
	N = nrow(X)
	D = ncol(X)
	Z = numeric(N)
	K = 0
	alpha = 1e-1

	itmax = 500

	# 3-3. execute gibbs sampling
	for(it in 1:itmax){

		# randomize samples
		idx = sample(1:N, N)

		# sample z for each x
		for(n in idx){

			Z[n] = 0
			x = t(X[n,])

			lnLn = numeric(K+1)
			if(K > 0){
				for(k in 1:K){
					dat_idx = which(Z == k)
					Nz = length(dat_idx)
					if(Nz > 0){
						lnLn[k] = log(Nz) + ln_marginal_likelihood(x, posterior(X[dat_idx,], theta0))
					}else{
						lnLn[k] = -Inf
					}
				}
			}
			lnLn[K+1] = log(alpha) + ln_marginal_likelihood(x, theta0)
			
			Z[n] = sample(1:(K+1), 1, prob=exp(lnLn - max(lnLn)))

			if(Z[n] == (K+1)){
				K = K + 1
			}

		}

	}	
}

無限混合正規分布

と書くと、とても格好良い。
最終兵器感がある。

別名、正規分布版のDirichlet Process Mixture(DPM)。ノンパラメトリックベイズ推定の1つである。

詳細は、M.D. Escober, M. West, “Bayesian Density Estimation and Inference Using Mixtures”やRadford M. Neal, “Markov Chain Sampling Methods for Dirichlet Process Mixture Models”に詳しい。

連続の確率分布を無限の点で離散化するDirichlet Processを混合モデルの事前分布に置くことで、理論的には無限の混合が可能になる・・・らしい。原理を理解するには、測度論という抽象数学が必要かつ自分にはハードルが高すぎなので、諦めてそういうものかと納得して使っている。

DPMは、実利的には、データの複雑さに応じてモデルの混合数を自動的に推定するための枠組みである。実装の詳細は、別の機会にまとめるが、Neal’s algorithm 3 (collapsed gibbs sampler) を正規分布の尤度関数+Normal-inverse-Wishartのパラメータ分布に対して、実行するとこんな感じになる。

iGM

左下の図が、分割数の学習による変化を表したグラフである。最初は単純な混合モデルに始まり、次第にデータに適した形で複雑さを学習していっていることがわかる。

通常の混合正規分布推定の場合、混合数は予め与えられていなければならない。混合数として、データ数Nを指定しておけば、同じことは一応できるが、全要素について計算しつづけなければならず、負荷は非常に大きい。一方、DPMの場合は、現時点の混合数K+1のみを扱えばよいので、遥かに効率的である。

個人的にはDeep Learningのようにブラックボックスなモデルよりも、こういうわかりやすい(?)モデルのほうが好き。

とある童話の統計科学(2)

Toaru_Logo

今回は、語の基本的な統計情報についての紹介。


語の統計

まずは各文書の語数とボキャブラリ数(ユニークな語の数)。

Toaru2_1

4文書合わせての総ボキャブラリ数は、11017語であった。

不思議の国のアリスは、完全に子供向けなのか(個人的にはあまりそうは思えないけれど。。。)、ボキャブラリ・語数共に控えめである。クリスマスキャロルは、語数は同程度だが、利用しているボキャブラリは1000語以上多い。low teen手前くらいだろうか。

次に、文書ごとの語の出現確率を示す。横軸は、各文書ごとに出現確率の高い順から語をソートし、その結果のランクを示している。縦軸は、各ランクの出現頻度のlog10である。文書ごとに、各ランクに対応する語は異なることに注意。

Toaru2_2

実はこのプロットから、「出現確率がランクの逆数に比例する」というZip則が概ね成立していることがわかる。適当にチューニングしてみると、出現確率=0.09/ランク、程度になる。どの文書でもあまり変わらない。

Zip則はロングテールの存在を表しており、ランクに対する累積確率を下図のように計算していくと、後ろのランクまで1.0に収束することがない。文書における語の生成は、基本的に疎と言える。

Toaru2_3

ちなみに、累積確率の最初の立ち上がりが非常に早いが、これらは主に情報の少ない機能語(冠詞, 前置詞, etc)が占めている。不思議の国のアリスの場合、語の頻出トップ10は以下である。

  • $(終端記号) 20.2%
  • the 4.8%
  • and 2.5%
  • to 2.1%
  • a 1.8%
  • it 1.7%
  • she 1.6%
  • i 1.6%
  • of 1.5%
  • said 1.3%

文書を特徴づける情報はロングテールに含まれており、連続データと同じ感覚で少数派をネグると痛い目をみる。


次回は、統計モデルn-gramのお話。

とある童話の統計科学(1)

Toaru_Logo

TeXでとある魔術のロゴ作成 というページを見つけたのがきっかけで、今回のネタをやってみようと思った(texソースお借りしました。ありがとうございます。素晴らしい出来ですね)。

何がしたいかというと、童話を対象に統計的自然言語処理をしてみようというお話。n-gram, LDA, HMMとフォローしていく予定だけど、長くなりそうなので、少しずつ。どこまでいけるかな。

今回はイントロダクションとして、何を対象にしたのかと文字の統計についての紹介。


準備

対象とする童話は、

  • Alice’s Adventures in Wonderland
  • A Christmas Carol
  • The Adventures of Tom Sawyer
  • Grimms’ Fairy Tales

いずれも、Project Gutenbergから入手した英語版。英語の方が、単語の切れ目がはっきりしていて扱いやすい。選択に特に意味はない。

前処理

今回はあまり関係ないけれど、主な前処理は以下。

  • Project Gutenbergのヘッダ・フッタの除去
  • 小文字化
  • [a-z0-9,.!?’]以外の文字を空白で置換
  • 数字列の置換 (e.g. 123 -> #)
  • ‘,.!?’を終端文字’$’として置換

文字(アルファベット)の統計

まず、各文書の総文字数を表したグラフを示す。

Toaru1_1

不思議の国のアリスが最も少なく、グリム童話が最も多い。その差は、300k文字程度である。

次に、各文字の出現確率(各文字の出現回数を総文字数で割ったもの)を表したグラフを示す。

Toaru1_2

自然言語処理を勉強したときに、一番最初に衝撃を受けたのは、文書が異なっても文字の頻度分布にそれほど大きな差がないということだった。文書のサイズが大きく違っているにもかからず、上図の結果はそれを良く表している。


今回は、ここまで。統計的自然言語処理は、このような言語に潜むパターンをモデル化して、利用する研究分野である。

一通り勉強はしたけれど、忘れかけているので、少しずつ更新していってリマインダにする予定。

Gaussian-Bernoulli Restricted Boltzmann Machine

先のCD learningの考察によって、ようやく意味がわかって、入力が連続値タイプのGaussian-Bernoulli RBMを実装できた。ちょこちょこ動いて、とってもかわいい!

GRBMDemo

Hinton先生のありがたいPracticeはほとんど反映されていない超簡易実装R版。混合分布推定に比べるとかなりシンプルで驚いた。

# Gaussian Bernoulli Restricted Boltzmann Machine

library(mvtnorm)

# 1. data generation
N = 200
D = 2

mu = list(NULL)
mu[[1]] = c(-10,-10)
mu[[2]] = c(-10,0)
mu[[3]] = c(-10,10)
mu[[4]] = c(0,-10)
mu[[5]] = c(0,0)
mu[[6]] = c(0,10)
mu[[7]] = c(10,-10)
mu[[8]] = c(10,0)
mu[[9]] = c(10,10)

S = diag(D)

K = length(mu)
X = matrix(0, N, D)
Z = numeric(N)

for(n in 1:N){
	Z[n] = sample(1:K, 1)
	X[n,] = rmvnorm(1, mu[[Z[n]]], S)
}

plot(X[,1], X[,2])

# 2. GRBM modeling

# 2-1. set up data, hyperparameters and function
v0 = scale(X)    # data set
L = 200          # number of hidden nodes
eta = 1e-2       # learning rate
itmax = 200      # max number of iteration

sigmoid = function(x){
	1.0/(1.0 + exp(-x))
}

# 2-2. set up parameters & variables
W = matrix(rnorm(D*L, 0, 1e-2), D, L)
b = rnorm(D,0,1e-2)
c = rnorm(D,0,1e-2)
R = diag(D) * 1e-2

# 2-3. execute learning
par(mfrow=c(2,2))
par(cex=0.6)
par(mar=rep(2,4))

rss = rep(NA, itmax)

for(it in 1:itmax){

	# a. calculate reconstructed samples and probability of hidden nodes; h0, v1, h1
	h0 = matrix(0, N, L)
	v1 = matrix(0, N, D)
	h1 = matrix(0, N, L)

	Rinv = solve(R)

	for(n in 1:N){
		h0[n,] = as.double(runif(L) < sigmoid(c + t(v0[n,]) %*% Rinv %*% W))
		v1[n,] = rmvnorm(1, b + W %*% h0[n,], R)
		h1[n,] = sigmoid(c + t(v1[n,]) %*% Rinv %*% W)	
	}
	
	# b. calculate statistics
	Evht0 = t(v0) %*% h0
	Evht1 = t(v1) %*% h1
	Ev0 = apply(v0, 2, mean)
	Ev1 = apply(v1, 2, mean)
	Eh0 = apply(h0, 2, mean)
	Eh1 = apply(h1, 2, mean)

	Evvt0 = Ev0 %*% t(Ev0) + cov(v0)
	Evvt1 = Ev1 %*% t(Ev1) + cov(v1)

	rss[it] = sum((v0 - v1)^2)

	# c. calculate derivatives of the contrastive divergence (CD)
	dW = - Rinv %*% (Evht0 - Evht1)
	db = - Rinv %*% (Ev0 - Ev1)
	dc = - (Eh0 - Eh1)

	tmp1 = -0.5 * Rinv %*% (Evvt0 - 2*Ev0%*%t(b) - 2*Evht0%*%t(W)) %*% Rinv
	tmp2 = -0.5 * Rinv %*% (Evvt1 - 2*Ev1%*%t(b) - 2*Evht1%*%t(W)) %*% Rinv
	dR = diag(diag(tmp1 - tmp2))

	# d. gradient-descent update
	W = W - (eta * max(abs(W)) / max(abs(dW))) * dW
	b = b - (eta * max(abs(b)) / max(abs(db))) * db
	c = c - (eta * max(abs(c)) / max(abs(dc))) * dc
	R = R - (eta * max(abs(R)) / max(abs(dR))) * dR

	# e. plot the progress
	plot(v0[,1], v0[,2], col="gray", main=sprintf("it=%d, v0=gray, v1=blue",it), xlim=c(-2,2), ylim=c(-2,2))
	points(v1[,1], v1[,2], col="blue")

	plot(1:itmax, rss, type="l", main=sprintf("it=%d, rss",it), ylim=c(0,200))
	abline(h=36, col="red")
}


rev.1 プログラムにミスがあったので修正。勾配法の更新量を調整して安定化。

MLE vs. Contrastive Divergence Learning

Restricted Boltzman Machineのプログラムを組んでみたのだけれど、上手く学習ができない。

真面目にCD learningの勉強をしなければと思い立ち、1次元ガウスの場合で計算をしてみた。以下、その結果。


下図は、ガウスの平均と精度パラメータの学習過程である。

KLvsCD_fig1

Gradient descent版の最尤推定と比較してみると、学習が遅いことがわかる。しかし、下図のように尤度でみると、毎イテレーションで、基本的に単調増加はしている。ただし、保証はされていない。

KLvsCD_fig2

おまけに、今回真面目に計算して分かったのだが、Contrastive Divergenceとは経験分布(データ)と真のモデル分布とのKullback Leibler Divergenceの超荒い近似であるようだ。RBMでは、潜在変数のGibbs Samplingが必要なため、さらにその近似精度が下がるであろうことは推測できる。

すなわち、学習には忍耐が必要・・・ということか。

上記の詳細はこちら


ということで、CD学習が正しく機能することは分かったので、もう少し頑張ってみようかな。

ノート on iPad

20150425a

iPadでノート取れるようになったかな、と久しぶりに思い立ち、
iOSのNoteshelfというアプリと
Jot Pro 2.0というスタイラスでトライしてみた。

2〜3年前に検討したときのiPadでノートを取ろう計画の満足度は10/100点。
まず、細い字がかけるようなスタイラスはなかった(と思う)し、
アプリの筆記入力に対するレスポンスも遅すぎた。

けど、今回は、及第の60点くらい。自分的には、実用できるレベルになった。
細い字も描きやすいし、入力レスポンスもかなり早くなった(若干ディレイがある)。

次の2〜3年は、レスポンスの高速化と、スタイラスの軽量化に期待!

Inverse-Wishart分布からInverse-Gamma分布への還元

Wishart分布がGamma分布の多変量拡張なので、Inverse-Wishart分布からInverse-Gamma分布への還元が成立するのは当然なのだが備忘録として、この変形を残しておく。

IW2IG

詳細はこちら

変形するとわかるが、Wishart分布とGamma分布では有効観測数(degree of freedom)の重みが2倍違っている。
つまり、同じ有効観測数のpriorを考えたときに、Inverse-Wishart分布のpriorとしての重みは、Inverse-Gamma分布のpriorとしての重みよりも小さい。


同じ有効観測数なのにと思って使っていたら、Wishart分布の挙動が不安定すぎて気づいた事実でした。