boost-python ではじめる大規模機械学習 （２） - データの受け渡し

あらすじ

boost-python を使用して、Python と C 言語両方を活用する方法を説明しています。
前記事では、C 言語で Hello world を出力する関数を作成し、これを Python から呼び出す方法を説明しました。
本記事では、ふたつの言語の間で簡単なデータを受け渡す方法について解説します。

プリミティブの受け渡しは簡単

渡したいデータが文字列や数値である場合は、何も気兼ねする必要がありません。

hello.cpp

#include <string>
#include <boost/python.hpp>
using namespace boost::python;

double increment(int x)
{
    return x + 1;
}

std::string add_txt(std::string x)
{
    return x + ".txt";
}

BOOST_PYTHON_MODULE(hello)
{
    def("increment", increment);
    def("add_txt", add_txt);
}

実行結果

実行時型チェックもあります。すばらしいです。

$ python
>>> import hello
>>> hello.increment(2)
3.0
>>> hello.add_txt('test')
'test.txt'
>>> hello.increment()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
Boost.Python.ArgumentError: Python argument types in
    hello.increment()
did not match C++ signature:
    increment(int)
>>>

タプルの受け渡しも簡単

Python のタプルも簡単に扱うことができます。
extract 関数を使用して、Python オブジェクトを C 言語のプリミティブに変換します。

hello.cpp

#include <string>
#include <boost/python.hpp>
using namespace boost::python;

object export_tuple()
{
    return make_tuple(1, 2, 3);
}

void import_tuple(object tup)
{
    size_t num = len(tup);

    for (size_t i = 0; i < num; i++)
    {
        std::cout << (extract<double>(tup[i]) + 1) << std::endl;
    }
}

BOOST_PYTHON_MODULE(hello)
{
    def("export_tuple", export_tuple);
    def("import_tuple", import_tuple);
}

実行結果

C 言語側では、tup オブジェクトにインデックスでアクセスしているだけであり、引数がタプルであるとは書いていません。したがって、リスト型でも動作するコードになっています。
型安全とはどこへやらという感じです。Python らしい。

>>> import hello
>>> hello.export_tuple()
(1, 2, 3)
>>> hello.import_tuple((1, 2, 3))
2
3
4
>>> hello.import_tuple([1, 2, 3])
2
3
4
>>> 

次回からが本番です。続く...

boost-python ではじめる大規模機械学習 （１）

はじめに

Python は、機械学習の分野で広く使用されるスクリプト言語です。SciPy や matplotlib といった、科学計算に特化したライブラリが多数提供されているのが特徴です。
いっぽう、弱点もあります。for ループの速度が遅いこと、並列処理が苦手なことなどです。これらの目的には、C 言語が適しています。
そこで、本記事では、Python と C の相互連携を可能にする boost-python ライブラリを使用して、大規模科学計算を効率的に解く方法を紹介します。題材には混合ガウス分布を使用します。
動作テストは、Fedora 環境で行なっています。Ubuntu でも動くと思います。

準備

はじめに、SciPy、matplotlib、boost ライブラリをインストールしてください。
boost-devel とか色々インストールした記憶がありますが、忘れました…

Hello, world!

最初はやっぱり Hello world

hello.cpp

#include <iostream>
#include <boost/python.hpp>
using namespace boost::python;

void hello()
{
    std::cout << "Hello, world!" << std::endl;
}

// 以下で定義した名前が、Python 側に公開されます。
BOOST_PYTHON_MODULE(hello)
{
    def("hello", hello);
}

setup.py

C++ で書いた拡張モジュールをコンパイルするためのコードです。
boost と OpenMP を読み込むよう設定しています。
GNU Scientific Library (GSL) なども読み込めます。

import sys
from distutils.core import setup, Extension

modules = []
modules.append(['hello', 'hello.cpp'])

extensions = []
for module in modules:
    mod_name = module[0]
    mod_source = module[1]
    ex_module = Extension(
        mod_name,
        sources = [mod_source],
        include_dirs = ['/usr/include/boost'],
        libraries = ['boost_python'],
        extra_compile_args=['-fopenmp', '-O3'],
        extra_link_args=['-lgomp', '-lpthread']
    )
    extensions.append(ex_module)

setup(name = "my-library", version = "1.0", ext_modules = extensions)

コンパイル

-i オプションで、同じディレクトリにコンパイルを行います。(付けないとよく分からないフォルダができる)
-f オプションは強制上書きです。

$ python setup.py build_ext -i -f
running build_ext
building 'hello' extension
...
gcc ...
g++ ...
$

これを毎回打つのが結構手間なので、

alias pyc='python setup.py build_ext -i -f'

としてます。
これで、作業ディレクトリに hello.so というバイナリコードが生成されます。

test.py

コンパイルしたモジュールは、Python モジュールと同様に読み込むことができます。

import hello
hello.hello()

実行結果

$ python test.py
Hello, world!
$

次回へ続く...

Roland Jupiter-8 解説

20世紀の音楽シーンは、Wurlitzer, Rhodes を始めとする様々な電子楽器により支えられてきました。その中でも、Roland が 1981 年に発売した Jupiter-8 は日本が産み出した最高傑作のひとつに数えるべきでしょう。本記事では、そんな Jupiter-8 についてゆるゆると解説したいと思います。

This photo is created by Ed Uthman, and licensed under Creative Commons (CC BY-SA 2.0).
http://www.flickr.com/photos/78147607@N00/2039658601/

時代背景…

大容量サンプラーが当たり前となった今日では想像もつかないことですが、1980 年頃の電子楽器最大の課題は同時発音数でした。ちなみに本機の同時発音数は 8 です。
ちなみに、この同時発音数の問題は、20 世紀の終わりごろまでついて回ることになります。
たとえば、1991 年発売、GM 規格の標準音源として有名な Roland SC-55mkII は 16 パート、同時発音数 28 でした。したがって、常に発音数を気にしておく必要がありました。

当時、発音数で世界を圧倒的にリードしていたのはヤマハでした。同時に複数の音が鳴らせるポリフォニック・シンセサイザーを実用化したのはヤマハでしたし、Jupiter-8 と同時代に発売された DX7 は同時発音数 16 と、両手で押さえきれない（！）数の音を同時に鳴らすことができました。

構造

Jupiter-8 はアナログシンセです。つまり、電子回路を物理的に発振させて音を作り、フィルター回路等で加工をしていきます。発振回路のことをオシレーターと呼びます。
パネル上のほぼ全てのスイッチの値はメモリー上に記録されており、64種類のプリセットを保存することが可能。これも当時としてはかなり画期的な仕掛けだった。
鍵盤は 61 鍵。

拡大写真: Roland Jupiter-8 Synth, 1983

上段

• ロゴ下
  • VOLUME
  • BALANCE (LOWER - UPPER)
  • ARPEGGIO
    • RATE
    • INT / EXT
• LFO
  • RATE, DELAY TIME, WAVE FORM のスイッチが見えます。大体想像の通りです。WAVE FORM ツマミで RANDOM が選択できるようです。
• VCO MODULATOR
  • （・ω・）…
• VCO-1
• VCO-2
  • VCO-1, 2 から VCA まで、細い横線が引かれています。この順で回路をとおる？
  • SOURCE MIX と書かれたダイアルが付いています。オシレータの信号を好きなバランスで混ぜられるということでしょうか。
• HPF
  • ハイパスフィルターを通るようです。カットオフ周波数しかありません。なんとも豪快です。
• VCF
  • フィルター回路は、VCF (Voltage-controlled filter) と VCA (Voltage-controlled amplifier) の二種類の回路で構成されています。VCF では、カットオフ周波数を操作します。
  • レゾナンスは設定できそうです。
• VCA
  • VCA では、音量を操作します。LFO 信号を入力して、わんわんと音量を上げ下げすることができます。
• ENV-1
• ENV-2
  • クラシカルな ADSR エンベロープが付いています。
  • 極性スイッチが付いているので、エンベロープを逆にできるようです。

下段

• MASTER TUNE
• TUNE
• ARPEGGIO
  • 1 / 2 / 3 / 4 / UP / DOWN / U&D / RANDOM
• ASSIGN MODE
  • SOLO / UNISON / POLY (1) / POLY (2)
• HOLD
  • LOWER / UPPER
• KEY MODE
  • DUAL / SPLIT / WHOLE
• PANEL MODE
  • LOWER / UPPER
• PATCH NUMBER
  • LOWER / UPPER
  • 画面中央に四桁のディスプレイがあります。鍵盤の上下で別々の音色が使えるようです。
• PATCH NUMBER 11-88
  • 1 / 2 / 3 / 4 / 5 / 6 / 7 / 8 / MANUAL
• PATCH PRESET
  • A / B / C / D / E / F / G / H / MEMORY PROTECT / WRITE
• TAPE MEMORY
  • SAVE / VERIFY / LOAD / DATA CHECK
  • このスイッチがなんとも時代を感じさせます。

鍵盤横

VCO とか LFO とかのスイッチがたくさんあります。

音色

甘く透き通るような音色が特徴的です。非常に日本的な音だとおもいます。
ハイがしなしなとしていて中低音がもやっとした出音には、ローランドの原点が垣間見えます。
それにしてもサインリード好きな会社ですね。

中段のみ Jupiter-8。他の鍵盤はジュピターではないので注意。

プリセットが色々聴けます。パッドやエレピ、ストリングの音など、これだけの回路の割にかなり色々な音が鳴ります。やはり高音は伸びません。

LFO 横に 2 つ LED ぽい何かがあります。これ光るのね…
それにしてもこの外人ノリノリである

参考

Roland Jupiter-8 - Wikipedia, the free encyclopedia
ローランド・SCシリーズ - Wikipedia
proun.net - ROLAND JUPITER-8

結局、すべては基本

早起きと適正体重に勝るものなし。なのかな。

五年後の自分、とは？

「五年後の自分について説明してください。」
エントリーシートで、こう問われることがあります。
一見回答不可能なこの問題は、実は「機会をどう選別するか」という問いなのだと思うようになりました。

まず、ストレートに考えます。
五年後の自分について、客観的に説明できるだろうか。
現在から因果の糸をたどって、五年後の自分にたどり着くだろうか。

無理だと思います。
五年前に想像した、五年後の未来では、私は理学研究科で修士二年で、素粒子物理学の研究をしているはずでした。
実際は、情報学研究科で修士一年で、音楽の統計解析をしています。
昔なら考えもしなかったことです。

五年後の未来を論理的に予測するのは、そもそもバカみたいな話です。
そういう、未来に対する論理的な予測能力は、経歴を見て面接をすれば、すぐに分かりますし、「五年後の自分」という最も非論理的な設問では測れません。

手掛かりは、「主体性」と「機会」なのではないか、と思うようになりました。
人は主体性によって機会を創り出し、機会によって成長します。主体性があるからこそ、結果に責任を持つことができます。
リクルートの（昔の）経営理念である、「自ら機会を創り出し、機会によって自らを変えよ」という言葉が参考になると思います。

やりたいこと、挑戦したいことを、どんどん試せばいい。
でもいつか時間的な限界にぶつかります。
そもそも、機会なんて、掃いて捨てるほどあるのです。
Twitter で誰かを誘ってご飯を食べたことがないなら、誘えばいいし、カラオケの幹事経験がなければ、やればいいのです。
何でも企画して、実行すれば良いのです。
そのうち、機会というものがどんどん見えてくるようになるし、全部やるのは無理、というのも分かります。

どういう機会を選ぶのか。
選択の基準はどこにあるのか。
それを、五年後、なりたい自分から逆算して考えろ、ということなのかも。

電子楽器の歴史

20 世紀の音楽を創り上げてきたあれやこれやの電子楽器が、あまり知られていないようなのでまとめてみました。
おまけとして、日本メーカー各社の最新モデルも載せました。
半世紀を通して、ハモンドやローズなど、物理的な発音システムの音をアンプで大きくするものから、真空管や IC で音を組み立てていくシンセサイザーへの移行が進んでいきます。アナログシンセについては、挙げ始めるとキリがないので、ちょっとだけ載せました。

※動画ごとに音量がバラバラなので、注意して聴いてください。
※それぞれの動画に関する権利は、各動画の著作者に帰属します。

Hammond Organ

1940 年代のロックシーンで爆発的人気を誇ったハモンドオルガンとレスリースピーカー。このスピーカーは、ツイーターとウーファーそれぞれが回転するようになっていて、鍵盤左手に取り付けられたハーフムーンスイッチを切り替えることで回転数を変化させることができる。

Rhodes Mark II

いわゆる、エレクトリックピアノ。金属製の音叉をハンマーでたたき、電磁場の変動を音として取り出す。Mark II は 1980 年頃のモデルで、音が強く歪んでいるのが特徴。上段に積まれているのは KORG SV-1。

Rhodes Mark V

これも有名なモデル。1984 年発売。メロウな音が特徴的。ただし、Mark II でもこれに近い音は鳴る。Mark V 発売当時、ダイノマイローズ（ダイノマイピアノ）と呼ばれるきらびやかな音色セッティングが流行ったが、この音が YAMAHA DX7 で再現しやすかったこともあり、ローズピアノの衰退を早めたと言われる。

Wurlitzer 140B

Rhodes と並び、エレクトリックピアノの音色のイメージを作り上げた。カーペンターズが好んで使用したことで有名。

Sequential Circuits Prophet 5

1978 年発売。代表的アナログ・シンセサイザーで、Yellow Magic Orchestra の曲などによく使われている。実は 5 音しか鳴らせない。当時はシンセサイザーの回路が高価で、同時に何音鳴らせるかというのが結構重要な課題だった。

Roland Jupiter-8

1980 年発売。アナログの電子回路を組み合わせて音を作っていく。ハモンドオルガンやローズピアノと異なり、物理的な振動装置が入っていないのが特徴。甘く染みこむような音色が特徴的。リズムを鳴らしているのは Roland TR-808。

YAMAHA DX7

1983 年発売。CPU などの電子回路を使って音を鳴らす、デジタル方式のシンセサイザーの原点として知られるモデル。

Roland SC-88pro

1996 年発売。安価なシンセサイザーの登場とインターネットの普及で、誰でも手軽に自分の音楽を公開できるようになった。

YAMAHA MOTIF XF

2010 年発売。YAMAHA 製の最新シンセサイザー。

Roland GR-55

2011 年発売の、ローランド製ギターシンセサイザー。ギターの音を入力にして、様々な音色を作ることができる。

KORG KAOSS PAD QUAD

2011 年発売。
おなかいっぱいです

おまけ：Wurlitzer SideMan

これ、リズムマシンなんだぜ…

動く変分混合ガウス分布（実装編）- 動く PRML シリーズ（２）

こちらもどうぞ - 動く変分混合ガウス分布（導出編）

実装には python, SciPy と matplotlib を使います。
テストデータには Old Faithful 間欠泉データを使います。
また、データの読み込み、プロットは混合ガウス分布の際に実装したものを再利用しますので、こちらからダウンロードしておいてください。

必要な関数の読み込み

はじめに、必要な関数を読み込みます。

from gmm import faithful_norm, init_figure, preview_stage
from scipy import arange, array, exp, eye, float64, log, maximum, ones, outer, pi, rand, zeros
from scipy.linalg import det, inv
from scipy.maxentropy import logsumexp
from scipy.special import digamma

主な関数は以下の通り。

• faithful_norm … データの読み込み
• init_figure … グラフウィンドウの作成
• preview_stage … 混合ガウス分布のレンタリング
• arange … 非負数列 [0, 1, 2, …] の作成
• eye … 単位行列の作成
• outer … テンソル積 の計算
• logsumexp … の計算
• digamma … ディガンマ関数

関数定義、データ読み込み

前回と同じく、estimate 関数を実行するとデモが見られるようにしておきます。

def estimate(num_class = 6, num_iter = 100, alpha0 = None, beta0 = None,
             m0 = None, W0 = None, nu0 = None):
    # Data
    x = faithful_norm()
    num, ndim = x.shape

    # Graph
    ax = init_figure()

今回は変分ベイズ法を使用するため、クラス数は多めに指定しておくと自動的に最適化されます。今回は 6 クラスを指定しています。繰り返し回数は 100 回としました。
クラス負担率、平均、精度のハイパーパラメータ は、関数呼び出しの際に自由に値を設定できるようにしておきます。明示されていない時 (None の時) は後でデフォルト値を代入します。

ハイパーパラメータの初期化

None の場合のみデフォルト値に差し替える関数を作ります。

def fill_param(param, default):
    if (param == None):
        return default
    else:
        return param

次に、ハイパーパラメータの値を設定します。デフォルト値は無情報事前分布 (non-informative prior) と呼ばれる値で、これは推論に恣意的な事前情報を持ち込まないことを意味しています。（alpha0 = 1000 とかで動かしてみるとよく分かります。）
ただし、精度に無情報事前分布を使うと不安定になりがちなので、精度 1，有効観測数 1 の事前情報を与えます。

# Hyperparameters
alpha0 = fill_param(alpha0, ones(num_class) * 1e-3)
beta0 = fill_param(beta0, 1e-3)
m0 = fill_param(m0, zeros(ndim))
W0 = fill_param(W0, eye(2))
nu0 = fill_param(nu0, 1)
inv_W0 = inv(W0)

負担率の初期化

今回は変分 E ステップで初期化を行うので、負担率を一様乱数からサンプリングします。

# Initial VB-E Step
r_nk = rand(num, num_class)
r_nk /= r_nk.sum(1)[:, None]

変分 M ステップ

負担率が手元にあるので、この値を使ってハイパーパラメータを更新していきます。

クラス混合比の更新

　　

alpha = alpha0 + Nk

平均、精度の更新

変分ベイズ法の更新式は、基本的に最尤推定の結果と事前分布の重み付き平均になるので、はじめに最尤推定の場合の結果を計算します。
　　

def calc_xbar(x, r_nk):
    num, ndim = x.shape
    num, num_class = r_nk.shape
    ret = zeros([num_class, ndim], dtype = float64)

    for k in xrange(num_class):
        clres = r_nk[:, k]
        for i in xrange(ndim):
            ret[k, i] = (clres * x[:, i]).sum()
        ret[k, :] /= clres.sum()

    return ret

　　
python の特性上、N でループを回すより K で回したほうが速いので、次の式を使って計算します。
　　

def calc_S(x, xbar, r_nk):
    num, ndim = x.shape
    num, num_class = r_nk.shape
    ret = zeros([num_class, ndim, ndim], dtype = float64)

    for k in xrange(num_class):
        clres = r_nk[:, k]
        for i in xrange(ndim):
            diff_i = x[:, i] - xbar[k, i]
            for j in xrange(ndim):
                diff_j = x[:, j] - xbar[k, j]
                ret[k, i, j] = (clres * diff_i * diff_j).sum()
        ret[k] /= clres.sum()

• 平均の更新

　　

def calc_m(xbar, Nk, m0, beta0, beta):
    num_class, ndim = xbar.shape
    ret = zeros([num_class, ndim], dtype = float64)

    for k in xrange(num_class):
        ret[k] = (beta0 * m0 + Nk[k] * xbar[k]) / beta[k]

    return ret

• 精度の更新

　　

def calc_W(xbar, Sk, Nk, m0, beta0, inv_W0):
    num_class, ndim = xbar.shape
    ret = zeros([num_class, ndim, ndim], dtype = float64)

    for k in xrange(num_class):
        ret[k] = inv_W0 + Nk[k] * Sk[k]
        fact = beta0 * Nk[k] / (beta0 + Nk[k])
        diff = xbar[k] - m0
        for i in xrange(ndim):
            for j in xrange(ndim):
                term = diff[i] * diff[j]
                ret[k, i, j] += fact * term
        ret[k] = inv(ret[k])

    return ret

最後に、ここまでの関数をつなぎ合わせて変分 M ステップを組みます。

for iiter in xrange(num_iter):
    # VB-M Step (Dirichlet)
    Nk = r_nk.sum(0)
    alpha = alpha0 + Nk

    # VB-M Step (Normal-Wishart)
    xbar = calc_xbar(x, r_nk)
    Sk = calc_S(x, xbar, r_nk)
    beta = beta0 + Nk
    mk = calc_m(xbar, Nk, m0, beta0, beta)
    W = calc_W(xbar, Sk, Nk, m0, beta0, inv_W0)
    nu = nu0 + Nk

    # VB-E Step
    ...
    # Visualization
    ...

変分 E ステップ

変分 E ステップの実装では、様々な確率変数の期待値を計算する必要があるので、これを順番に実装していきます。プロットに必要な期待値の計算も行います。モクモクと作ります。（PRML を参照のこと。）

ディリクレ分布

　　

def expect_pi(alpha):
    return alpha / alpha.sum()

　　

def expect_lpi(alpha):
    return digamma(alpha) - digamma(alpha.sum())

ウィシャート分布

　　

def expect_llambda(W, nu):
    ndim = W.shape[0]
    arr = float64(nu - arange(ndim)) / 2
    return digamma(arr).sum() + ndim * ln2 + log(det(W))

　　

def expect_lambda(W, nu):
    return W * nu

正規分布

　　
ここでは、二次形式 を計算する関数を実装し、これを活用します。

def quad(A, x):
    num, ndim = x.shape
    ret = zeros(num, dtype = float64)

    for i in xrange(ndim):
        for j in xrange(ndim):
            ret += A[i, j] * x[:, i] * x[:, j]

    return ret

def expect_quad(x, m, beta, W, nu):
    ndim = x.shape[1]
    return ndim / beta + nu * quad(W, x - m[None, :])

　　

def expect_log(x, m, beta, W, nu):
    ndim = x.shape[1]

    ex_llambda = expect_llambda(W, nu)
    ex_quad = expect_quad(x, m, beta, W, nu)
    return (ex_llambda - ndim * ln2pi - ex_quad) / 2

負担率の更新

導出編に示した通り、対数負担率は
　　
という形で更新できるので、これを実装します。

ex_lpi = expect_lpi(alpha)
ex_log = zeros([num, num_class], dtype = float64)
for k in xrange(num_class):
    ex_log[:, k] = expect_log(x, mk[k], beta[k], W[k], nu[k])
lrho = ex_lpi[None, :] + ex_log

さらに、データ点ごとの負担率の総和が 1 になるように正規化します。

def normalize_response(lrho):
    num, num_class = lrho.shape
    ret = zeros([num, num_class], dtype = float64)

    for i in xrange(num):
        lr = lrho[i] - logsumexp(lrho[i])
        ret[i] = exp(lr)

    ret = maximum(ret, 1e-10) # ゼロよけ
    ret /= ret.sum(1)[:, None]
    return ret

可視化

変分 E ステップ、変分 M ステップが組みあがったので、これをレンタリングします。
ただし、最尤推定や MAP 推定と違い、今手元には や の値がありません。そこで、これらの値の期待値を計算してなんとなくそれっぽいグラフを作ることにします。

ex_pi = expect_pi(alpha)
ex_lambda = zeros([num_class, ndim, ndim], dtype = float64)
for k in xrange(num_class):
    ex_lambda[k] = expect_lambda(W[k], nu[k])
preview_stage(ax, x, ex_pi, mk, inverse_matrices(ex_lambda))

ソースコード

できあがり！
なんとなくこんな感じで動かすとホクホクします

$ python
Python 2.7.1 (r271:86832, Apr 12 2011, 16:16:18) 
[GCC 4.6.0 20110331 (Red Hat 4.6.0-2)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import vbgmm
>>> vbgmm.estimate()

# -*- coding: utf-8 -*-
from gmm import faithful_norm, init_figure, preview_stage
from scipy import arange, array, exp, eye, float64, log, maximum, ones, outer, pi, rand, zeros
from scipy.linalg import det, inv
from scipy.maxentropy import logsumexp
from scipy.special import digamma

ln2 = log(2)
ln2pi = log(2 * pi)

# ========== ========== ========== ==========
#
#   Helpers
#
# ========== ========== ========== ==========
#
#   パラメータが null の場合，デフォルト値を返します．
#
def fill_param(param, default):
    if (param == None):
        return default
    else:
        return param

#
#   行列のリストから，対応する逆行列のリストを計算します．
#
def inverse_matrices(A):
    return array(map(inv, A))

# ========== ========== ========== ==========
#
#   Dirichlet Distribution
#
# ========== ========== ========== ==========
#
#   E[\pi_k] を計算します．
#
def expect_pi(alpha):
    return alpha / alpha.sum()

#
#   E[\ln \pi_k] を計算します．
#
def expect_lpi(alpha):
    return digamma(alpha) - digamma(alpha.sum())

# ========== ========== ========== ==========
#
#   Wishart Distribution
#
# ========== ========== ========== ==========
#
#   E[\ln \Lambda] を計算します．
#
def expect_llambda(W, nu):
    ndim = W.shape[0]
    arr = float64(nu - arange(ndim)) / 2
    return digamma(arr).sum() + ndim * ln2 + log(det(W))

#
#   E[\Lambda] を計算します．
#
def expect_lambda(W, nu):
    return W * nu

# ========== ========== ========== ==========
#
#   Normal Distribution
#
# ========== ========== ========== ==========
#
#   行列 A とベクトルの配列 {x1, x2, …, xN} から
#   二次形式の配列 {A[x1], A[x2], …, A[xN]} を計算します．
#
def quad(A, x):
    num, ndim = x.shape
    ret = zeros(num, dtype = float64)

    for i in xrange(ndim):
        for j in xrange(ndim):
            ret += A[i, j] * x[:, i] * x[:, j]

    return ret

#
#   E[\lambda (x - \mu) ** 2] を計算します．
#
def expect_quad(x, m, beta, W, nu):
    ndim = x.shape[1]
    return ndim / beta + nu * quad(W, x - m[None, :])

#
#   E[\ln N(x|\pi, \mu, \Lambda)] を計算します．
#
def expect_log(x, m, beta, W, nu):
    ndim = x.shape[1]

    ex_llambda = expect_llambda(W, nu)
    ex_quad = expect_quad(x, m, beta, W, nu)
    return (ex_llambda - ndim * ln2pi - ex_quad) / 2

# ========== ========== ========== ==========
#
#   VB-E Step Helper
#
# ========== ========== ========== ==========
#
#   対数負担率を正規化します．
#
def normalize_response(lrho):
    num, num_class = lrho.shape
    ret = zeros([num, num_class], dtype = float64)

    for i in xrange(num):
        lr = lrho[i] - logsumexp(lrho[i])
        ret[i] = exp(lr)

    ret = maximum(ret, 1e-10)
    ret /= ret.sum(1)[:, None]
    return ret

# ========== ========== ========== ==========
#
#   VB-M Step Helper
#
# ========== ========== ========== ==========
#
#   負担率を用いてクラス平均を計算します．
#
def calc_xbar(x, r_nk):
    num, ndim = x.shape
    num, num_class = r_nk.shape
    ret = zeros([num_class, ndim], dtype = float64)

    for k in xrange(num_class):
        clres = r_nk[:, k]
        for i in xrange(ndim):
            ret[k, i] = (clres * x[:, i]).sum()
        ret[k, :] /= clres.sum()

    return ret

#
#   負担率を用いてクラス分散を計算します．
#
def calc_S(x, xbar, r_nk):
    num, ndim = x.shape
    num, num_class = r_nk.shape
    ret = zeros([num_class, ndim, ndim], dtype = float64)

    for k in xrange(num_class):
        clres = r_nk[:, k]
        for i in xrange(ndim):
            diff_i = x[:, i] - xbar[k, i]
            for j in xrange(ndim):
                diff_j = x[:, j] - xbar[k, j]
                ret[k, i, j] = (clres * diff_i * diff_j).sum()
        ret[k] /= clres.sum()

    return ret


#
#   平均のハイパーパラメータを計算します．
#
def calc_m(xbar, Nk, m0, beta0, beta):
    num_class, ndim = xbar.shape
    ret = zeros([num_class, ndim], dtype = float64)

    for k in xrange(num_class):
        ret[k] = (beta0 * m0 + Nk[k] * xbar[k]) / beta[k]

    return ret

#
#   精度のハイパーパラメータを計算します．
#
def calc_W(xbar, Sk, Nk, m0, beta0, inv_W0):
    num_class, ndim = xbar.shape
    ret = zeros([num_class, ndim, ndim], dtype = float64)

    for k in xrange(num_class):
        ret[k] = inv_W0 + Nk[k] * Sk[k]
        fact = beta0 * Nk[k] / (beta0 + Nk[k])
        diff = xbar[k] - m0
        for i in xrange(ndim):
            for j in xrange(ndim):
                term = diff[i] * diff[j]
                ret[k, i, j] += fact * term
        ret[k] = inv(ret[k])

    return ret

# ========== ========== ========== ==========
#
#   Program
#
# ========== ========== ========== ==========
def estimate(num_class = 6, num_iter = 100, alpha0 = None, beta0 = None,
             m0 = None, W0 = None, nu0 = None):
    # Data
    x = faithful_norm()
    num, ndim = x.shape

    # Graph
    ax = init_figure()

    # Hyperparameters
    alpha0 = fill_param(alpha0, ones(num_class) * 1e-3)
    beta0 = fill_param(beta0, 1e-3)
    m0 = fill_param(m0, zeros(ndim))
    W0 = fill_param(W0, eye(2))
    nu0 = fill_param(nu0, 1)
    inv_W0 = inv(W0)

    # Initial VB-E Step
    r_nk = rand(num, num_class)
    r_nk /= r_nk.sum(1)[:, None]

    # Execution
    for iiter in xrange(num_iter):
        # VB-M Step (Dirichlet)
        Nk = r_nk.sum(0)
        alpha = alpha0 + Nk

        # VB-M Step (Normal-Wishart)
        xbar = calc_xbar(x, r_nk)
        Sk = calc_S(x, xbar, r_nk)
        beta = beta0 + Nk
        mk = calc_m(xbar, Nk, m0, beta0, beta)
        W = calc_W(xbar, Sk, Nk, m0, beta0, inv_W0)
        nu = nu0 + Nk

        # VB-E Step
        ex_lpi = expect_lpi(alpha)
        ex_log = zeros([num, num_class], dtype = float64)
        for k in xrange(num_class):
            ex_log[:, k] = expect_log(x, mk[k], beta[k], W[k], nu[k])
        lrho = ex_lpi[None, :] + ex_log
        r_nk = normalize_response(lrho)

        # Visualization
        ex_pi = expect_pi(alpha)
        ex_lambda = zeros([num_class, ndim, ndim], dtype = float64)
        for k in xrange(num_class):
            ex_lambda[k] = expect_lambda(W[k], nu[k])
        preview_stage(ax, x, ex_pi, mk, inverse_matrices(ex_lambda))

if (__name__ == '__main__'):
    estimate()