Loading web-font TeX/Math/Italic

December 31, 2009

PythonでLDAを実装してみる

Latent Dirichlet Allocationはテキストのような不連続データのための生成的確率モデル。入力はドキュメント。出力はドキュメントを特徴づける何か(tf-idfみたいなもん)。

基本的なアイディアは、あるドキュメントは潜在的ないくつかのトピックが混合していて、それぞれのトピックは語の分布で特徴づけられている、ということ。

論文[1]ではαとβというパラメータを用いてドキュメントが以下のように生成されると仮定している。
  1. ドキュメントのトピックの分布θがディリクレ分布Dir(α)に基づいて選ばれる。
  2. ドキュメントの語数N個になるまで以下を繰り返す。
    1. トピックznが多項分布Mult(θ)に基づいて選ばれる。
    2. 単語wnが確率p(wn|zn,β)で選ばれる。
ただし、トピックzの数をk個、単語wの種類をV個とすると、パラメータαはk次元のベクトル、βはk x V次元の行列でβij=p(wj|zi)。


Fig.1 LDAのグラフィカルモデル。Mはコーパス内のドキュメントの数。Nは各ドキュメントの語の数。αによって分布θが決まり、分布θに従ってzが選ばれ、zとβに従ってwが選ばれる。


ここで、ドキュメントにおけるαとβの値がわかれば、トピックがどんな割合であって(α)、そのトピックに関する語がどんな割合で存在するか(β)がわかる。つまり、ドキュメントが上のようなプロセスで生成されているとしてαとβの値はいくつかということを推定するのがLDAの目的。

αとβを推定する方法は変分ベイズEMアルゴリズムを利用するものやGibbs Samplerを利用するものなどが提案されています。また、いくつもの派生的なモデルも提案されています。

本稿では、論文[1]と
lda, a Latent Dirichlet Allocation package
http://chasen.org/~daiti-m/dist/lda/
のmatlabの方をパクっ参考にして、変分ベイズEMをPythonで実装してみました。

実行方法

要Numpy、SciPy。
python lda.py [-Nclasses] [-Iemmax] [-Ddemmax] [-Eepsilon] 入力ファイル名 出力ファイル名
classesはクラスの数。emmaxはEMの最大反復回数。demmaxは一つのドキュメントのEM最大反復回数。epsilonは収束条件。
たとえば、トピックの数が10個、「train」が入力で「model」に出力する場合は、
python lda.py -N10 train model

入力

テキストファイル。
1:1 2:4 5:2
1:2 3:3 5:1 6:1 7:1
2:4 5:1 7:1
ってな具合のフォーマット。行がドキュメント、<語のID>:<カウント>。SVMのライブラリでつかうフォーマットと似ている。上記の参考にしたサイトの上から1/3あたりのDownloadからcかmatlabバージョンをダウンロードして解凍するとtrainというファイルがありそのまま使える。

出力

<出力ファイル名>.alpha
<出力ファイル名>.beta
の二つ。alphaは長さがトピック数のリスト。betaはトピックの数×語の数の二次元リストになっている。はず。

References

[1] Blei et al., Latent Dirichlet Allocation, The Journal of Machine Learning Research, 2003.

ソース

lda.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import sys,getopt
from numpy import array,matrix,diag
from scipy import sum,log,exp,mean,dot,ones,zeros
from scipy.special import polygamma
from scipy.linalg import norm
from random import random
def usage():
print "usage: python lda.py [-Nclasses] [-Iemmax] [-Ddemmax] [-Eepsilon] train model"
sys.exit(0)
def main():
# Get params
try:
opts, args = getopt.getopt(sys.argv[1:], "N:I:D:E:h")
train = sys.argv[len(opts) + 1]
model = sys.argv[len(opts) + 2]
except:
usage()
sys.exit(0)
# Set params
k = 10
emmax = 100
demmax = 20
epsilon = 0.0001
for o, a in opts:
if o == '-N': k = int(a)
if o == '-I': emmax = int(a)
if o == '-D': demmax = int(a)
if o == '-E': epsilon = float(a)
if o == '-h': usage()
# Train
print train
train = open(train,'r').read()
alpha,beta = ldamain(train, k, emmax, demmax, epsilon)
# Write
print 'writing data..'
# alpha
print 'alpha -> %s.alpha' % model
writer = open(model + '.alpha','w')
writer.write(str(alpha.tolist()))
writer.close()
# beta
print 'beta -> %s.beta' % model
writer = open(model + '.beta','w')
writer.write(str(beta.tolist()))
writer.close()
def ldamain(train, k, emmax=100, demmax=20, epsilon=1.0e-4):
'''
Latent Dirichlet Allocation, standard model.
[alpha,beta] = ldamain(train,k,[emmax,demmax])
train : group of data ASCII text.
Each line represents a document and consists of pairs of <feature_id>:<count>.
e.g.
1:1 2:4 5:2
1:2 3:3 5:1 6:1 7:1
2:4 5:1 7:1
k : # of classes to assume
emmax : # of maximum VB-EM iteration (default 100)
demmax : # of maximum VB-EM iteration for a document (default 20)
epsilon: A threshold to determine the whole convergence of the estimation (default 0.0001)
'''
d = [ zip(*[ [int(x) for x in w.split(':')] for w in L.split()]) for L in train.split('\n') if L ]
return lda.train(d,k,emmax,demmax)
class lda():
'''
Latent Dirichlet Allocation, standard model.
[alpha,beta] = lda.train(d,k,[emmax,demmax])
d : data of documents
k : # of classes to assume
emmax : # of maximum VB-EM iteration (default 100)
demmax : # of maximum VB-EM iteration for a document (default 20)
'''
@staticmethod
def train(d, k, emmax=100, demmax=20, epsilon=1.0e-4):
'''
Latent Dirichlet Allocation, standard model.
[alpha,beta] = lda.train(d,k,[emmax,demmax])
d : data of documents
k : # of classes to assume
emmax : # of maximum VB-EM iteration (default 100)
demmax : # of maximum VB-EM iteration for a document (default 20)
'''
# # of documents
n = len(d)
# # of words
L = max(map(lambda x: max(x[0]), d)) + 1
# initialize
beta = ones((L, k)) / L
alpha = lda.normalize(sorted([random() for i in range(k)], reverse=True))
gammas = zeros((n, k))
lik = 0
plik = lik
print 'number of documents (n) = %d' % n
print 'number of words (l) = %d' % L
print 'number of latent classes (k) = %d' % k
for j in range(emmax):
print 'iteration %d/%d..\t' % (j+1, emmax),
#vb-esstep
betas = zeros((L, k))
for i in range(n):
gamma,q = lda.vbem(d[i], beta, alpha, demmax)
gammas[i,:] = gamma
betas = lda.accum_beta(betas,q,d[i])
#vb-mstep
alpha = lda.newton_alpha(gammas)
beta = lda.mnormalize(betas,0)
#converge?
lik = lda.lda_lik(d, beta, gammas)
print 'log-likelihoood =', lik
if j > 1 and abs((lik - plik) / lik) < epsilon:
if j < 5:
print
return lda.train(d, k, emmax, demmax) # try again
return
print 'converged'
return alpha, beta
plik = lik
@staticmethod
def vbem(d, beta, alpha0, emmax=20):
'''
[alpha,q] = vbem(d,beta,alpha0,[emmax])
calculates a document and words posterior for a document d.
alpha : Dirichlet posterior for a document d
q : (L * K) matrix of word posterior over latent classes
d : document data
beta :
alpha0 : Dirichlet prior of alpha
emmax : maximum # of VB-EM iteration.
'''
digamma = lambda x: polygamma(0,x)
L = len(d[0])
k = len(alpha0)
q = zeros((L, k))
nt = ones((1, k)) * L / k
pnt = nt
for j in range(emmax):
#vb-estep
q = lda.mnormalize(matrix(beta[d[0],:]) * matrix(diag(exp(digamma(alpha0 + nt))[0])), 1)
#vb-mstep
nt = matrix(d[1]) * q
#converge?
if j > 1 and lda.converged(nt, pnt, 1.0e-2):
break
pnt = nt.copy()
alpha = alpha0 + nt
return alpha, q
@staticmethod
def accum_beta(betas, q, t):
'''
betas = accum_beta(betas,q,t)
accumulates word posteriors to latent classes.
betas : (V * K) matrix of summand
q : (L * K) matrix of word posteriors
t : document of struct array
'''
betas[t[0],:] += matrix(diag(t[1])) * q
return betas
@staticmethod
def lda_lik(d, beta, gammas):
'''
lik = lda_lik(d, beta, gammas)
returns the likelihood of d, given LDA model of (beta, gammas).
'''
egamma = matrix(lda.mnormalize(gammas, 1))
lik = 0
n = len(d)
for i in range(n):
t = d[i]
lik += (matrix(t[1]) * log(matrix(beta[t[0],:]) * egamma[i,:].T))[0,0]
return lik
@staticmethod
def newton_alpha(gammas,maxiter=20,ini_alpha=[]):
'''
alpha = newton_alpha (gammas,[maxiter])
Newton-Raphson iteration of LDA Dirichlet prior.
gammas : matrix of Dirichlet posteriors (M * k)
maxiter : # of maximum iteration of Newton-Raphson
'''
digamma = lambda x: polygamma(0, x)
trigamma = lambda x: polygamma(1, x)
M,K = gammas.shape
if not M > 1:
return gammas[1,:]
if not len(ini_alpha) > 0:
ini_alpha = mean(gammas, 0) / K
L = 0
g = zeros((1,K))
pg = sum(digamma(gammas), 0) - sum(digamma(sum(gammas, 1)))
alpha = ini_alpha.copy()
palpha = zeros((1, K))
for t in range(maxiter):
L += 1
alpha0 = sum(alpha)
g = M * (digamma(alpha0) - digamma(alpha)) + pg
h = -1.0 / trigamma(alpha)
hgz = dot(h,g) / (1.0 / trigamma(alpha0) + sum(h))
for i in range(K):
alpha[i] = alpha[i] - h[i] * (g[i] - hgz) / M
if alpha[i] < 0:
return lda.newton_alpha(gammas, maxiter, ini_alpha / 10.0)
#converge?
if L > 1 and lda.converged(alpha, palpha, 1.0e-4):
break
palpha = alpha.copy()
return alpha
@staticmethod
def normalize(v):
return v / sum(v)
@staticmethod
def mnormalize(m, d=0):
'''
x = mnormalize(m, d)
normalizes a 2-D matrix m along the dimension d.
m : matrix
d : dimension to normalize (default 0)
'''
m = array(m)
v = sum(m, d)
if d == 0:
return m * matrix(diag(1.0 / v))
else:
return matrix(diag(1.0 / v)) * m
@staticmethod
def converged(u, udash, threshold=1.0e-3):
'''
converged(u,udash,threshold)
Returns 1 if u and udash are not different by the ratio threshold
'''
return norm(u - udash) / norm(u) < threshold
if __name__ == '__main__':
main()
view raw lda.py hosted with ❤ by GitHub


雑感

すごく…遅いです…

December 28, 2009

Pythonで順列や組み合わせを手に入れる

from itertools import ...

IteratorArguments
product()p, q, ... [repeat=1]
permutations()p[, r]
combinations()p, r

product('ABCD', repeat=2)AA AB AC AD BA BB BC BD CA CB CC CD DA DB DC DD
permutations('ABCD', 2)AB AC AD BA BC BD CA CB CD DA DB DC
combinations('ABCD', 2)AB AC AD BC BD CD

10.7. itertools — Functions creating iterators for efficient looping — Python v2.6.4 documentation

December 26, 2009

Chironで実行可能なSilverlightアプリをつくるときのコマンド

構成が
app/app.py
app/app.xaml
index.html
で、index.html内で
<param name="source" value="app.xap"/>
としたら
Chiron /d:app /z:app.xap

Google App Engineでテンプレートを使おうとしたらUnicodeDecodeError

Google App Engineでテンプレートを使おうとしたらエラーが出る。

UnicodeDecodeError: 'ascii' codec can't decode byte 0xe3 in position 116: ordinal not in range(128)

どうやらテンプレートで日本語を使おうとしたときのunicodeのエンコード・デコードのエラーっぽい。

…
sentence = 'すもももももももものうち。'
template_values = {'sentence':sentence}
path = os.path.join(os.path.dirname(__file__), 'index.html')
self.response.out.write(template.render(path,template_values))
…

これのtemplate.render()に

.decode('utf-8')

をつけて

…
sentence = 'すもももももももものうち。'
template_values = {'sentence':sentence}
path = os.path.join(os.path.dirname(__file__), 'index.html')
self.response.out.write(template.render(path,template_values).decode('utf-8'))
…

となおしたら、大丈夫だった。テンプレートはutf-8で。

December 25, 2009

MacでLaTeXはDrag & Drop pTeXがいい

Snow Leopardになっても小川版pTeXがすごくいい。

Drag & Drop pTeX - JIS X0212 for pTeX
http://www2.kumagaku.ac.jp/teacher/herogw/

ghostscriptのインストールも忘れずに!

  1. 古いものを削除。アンインストールツールが付いている。


  2. pTeX.appをアプリケーションフォルダにドラッグアンドドロップ。


  3. TeXShopのPath設定で(pdf)TeXを

    /Applications/pTeX.app/teTeX/bin

    に。TeX + dvips + distllerの設定も

    Xtexshop

    に。dvipdfmxでフォントを埋め込む。なお埋め込まない場合はXtexshop-ryu。文字コードの認識しない場合は

    Shift_JIS/ヒラギノフォント埋め込みの場合は“dotexshop”、同/フォント非埋め込みの場合は“dotexshop-ryu”
    EUC/ヒラギノフォント埋め込みの場合は“dotexshop-euc”、同/フォント非埋め込みの場合は“dotexshop-euc-ryu”
    UTF-8/ヒラギノフォント埋め込みの場合は“dotexshop-utf8”、同/フォント非埋め込みの場合は“dotexshop-utf8-ryu”


  4. ついでに文字コードもUTF-8に。


  5. あとは Typesetting の Default Script を Tex + DVI にして、Misc の BibTeX を jbibtexにする。

December 23, 2009

MacのEclipseでGoogle App EngineのPYTHONPATH




Eclipseで新しいGoogle App Engineのプロジェクトを作るときのGoogle App Engine Directory、あるいはプロジェクトを指定してからプロパティでPyDev - PYTHONPATHのExternal Librariesの場所は
/Applications/GoogleAppEngineLauncher.app/Contents/Resources/GoogleAppEngine-default.bundle/Contents/Resources/google_appengine
PyDev - PYTHONPATHのExternal LibrariesではAdd Source FolderしてからCtrl+Shift+Gで直接場所を指定できる。


/usr/local/bin/dev_appserver.py
"${project_loc}/src"
--port==9999

December 18, 2009

ことえりにZコマンドを



Google日本語入力で使える「Zコマンド」がすごいのようにことえりのユーザ辞書に単語を登録した。

z+[hjkl] = ←↓↑→
z+- = ~
z+[ = 『
z+] = 』
z+, = ‥
z+. = …
z+/ = ・

記号入のほうは登録するときは注意が必要。たとえば「z.」は「z。」で登録しないとうまく行かない。

ユーザ辞書にRegister Wordsするときのショートカットは

Ctrl + Shift + N

MacからEPDことEnthought Python Distributionを消す方法

以下引用。

• On Mac, you'll need to uninstall manually as Apple does not yet support a standard uninstall mechanism for .mpkg installers. The recommended uninstall commands are as follows. The instructions suffixed by a #* are specific to restoring a different Python interpreter, in this case MacPython 2.5 from python.org. If you want to restore a different Python, you'll need to know where your backup Python environment is and customize these commands based on that information.

cd /Library/Frameworks/Python.framework/Versions
sudo rm -rf 2.5.2001 Current
sudo ln -fhs 2.5 Current
cd /Applications
sudo rm -rf "EPD 2.5.2001"
cd /usr/local/bin
PYTHON=/Library/Frameworks/Python.framework/Versions/Current
sudo ln -fsh PYTHON/bin/python sudo ln -fsh PYTHON/bin/pythonw
sudo ln -fsh PYTHON/bin/python2.5 sudo ln -fsh PYTHON/bin/pythonw2.5

As an example, on OS X 10.6, to restore the default Python installation you would set the PYTHON variable as:

PYTHON=/System/Library/Frameworks/Python.framework/Versions/2.6
EPD - Frequently Asked Questions (FAQ) :: Products :: Enthought, Inc.

December 16, 2009

Ubuntu 9.10とWindowsのデュアルブートでデフォルトのOSをWindowsにする

UbuntuとWindows 7とのデュアルブートにしようとした。デフォルトで起動するOSをWindows 7にしたくてググったら、/boot/grub/menu.lstを書き換えましょう、という記事がたくさん見つかったのだけど、そもそもmenu.lstが見つからない。

で、さらに調べたらGrub2に関する記事を読め、と。

デフォルトのOSを切り替えるためにはStartUpManagerを使いなさいとのこと。

手順にしたがってソフトウェアセンターからStartUpManagerをインストールして、



システム > システム管理からStartUp-Managerを起動して設定した。


December 13, 2009

jQuery + PythonのJSONでのデータのやり取り


jQueryとPython CGIのデータのやり取りをJSONでおこないたい。
構成は、
index.html
sample.js
sample.py
(lib/jquery/jquery-1.3.2.js)

index.html

<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
"http://www.w3.org/TR/html4/loose.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
    <head>
        <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
        <title>jQuery + Python</title>
    </head>
    <body>
        <h1>jQuery + Python</h1>
    
  <form id="hoge">
   <input type="text" name="foo" />
   <input type="text" name="bar" />
   <input type="text" name="baz" />
  </form>
  
  <a id="foo" href="#">show</a>
  <a id="bar" href="#">clear</a>
  <div id="baz"></div>

  <script type="text/javascript" src="lib/jquery/jquery-1.3.2.js"></script>
  <script type="text/javascript" src="sample.js"></script>
    </body>
</html>


sample.js

p2jはパラメータからJSON形式に変換する関数。どこかのページに掲載されていたのですが、失念してしまいました。すいません。
var $j = jQuery.noConflict();

function p2j(d) {
 if (d.constructor != Array) {
  return d;
 }
 var data={};
 for(var i=0;i<d.length;i++) {
  if (typeof data[d[i].name] != 'undefined') {
   if (data[d[i].name].constructor!= Array) {
    data[d[i].name]=[data[d[i].name],d[i].value];
   } else {
    data[d[i].name].push(d[i].value);
   }
  } else {
   data[d[i].name]=d[i].value;
  }
 }
 return data;
};


$j(document).ready(function(){ 
 $j('#foo')
 .click(function() {  
  // クエリ
  var query = $j(":input").serializeArray();
  console.log(p2j(query));
  
  // GETリクエスト
  $j.get('sample.py', query, function(text) {
   // 結果の処理
   var json = JSON.parse(text);
   var html = "";
   for (var i in json) {
    html += json[i].key + ':' + json[i].value + ' ';
   }
   $j('#baz').html(html);   
  });
  
  return false;
 });
 
 $j('#bar')
 .click(function (){
  $j('#baz').html('');
 }); 
});


sample.py

JSON形式にするにはjsonモジュールかsimplejsonモジュールがいいかもしれません。
#!/usr/bin/env python
# -*- coding: utf-8 -*-

import cgi
import cgitb; cgitb.enable()
import json

print "Content-type: text/javascript; charset=utf-8"
print

form = cgi.FieldStorage()

foo = form.getfirst("foo", "")
bar = form.getfirst("bar", "")
baz = form.getfirst("baz", "")

print json.dumps([{'key':'foo','value':foo},
                  {'key':'bar','value':bar},
                  {'key':'baz','value':baz}])