二次元のガウス分布

去年の10月頃からPRML(Pattern Recognition and Machine Learning)を読んでいる。なかなか難しくて読み進めるだけでも骨が折れる。実務に使えるようになるにはまだ時間がかかるかもしれない（部分的に役に立っている部分はあるけれど）。

どこが難しいのか。色々あるけど、多次元空間への感覚のなさがネックになっていることが多いように思う。きっと結果がイメージできるためだろう、議論が一次元空間に閉じているうちは素直に頭に入ってくる。多次元空間を扱うことになった瞬間に「納得感」がなくなる。

この状況をなんとか打開したい。

いきなり一般の多次元空間に行くのは難しすぎるので、まず二次元空間に慣れることにする。特に、良くあらわれる二次元のガウス分布に慣れたい。そこで、pythonのmatplotlibを使って二次元のガウス分布をお絵描きして遊んでみた。

from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
import numpy.linalg as la
import functools
import sys

def calc_z(x, y, mu, cov):
    """calculates 2-dimensional gaussian"""
    bx = np.array([[x],[y]])

    coef = 1.0 / (2 * np.pi) / (la.det(cov) ** (0.5))

    exparg1 = -1.0 / 2.0 * np.transpose(bx - mu)
    exparg2 = la.inv(cov)
    exparg3 = bx - mu

    # use np.dot for matrix multiplication
    exparg = np.dot(np.dot(exparg1, exparg2), exparg3)

    return coef * np.exp(exparg)[0][0]

def gaussian3d(mu, cov):
    """plots 2-dimensional gaussian"""
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    xwidth = np.sqrt(abs(cov[0][0]))
    ywidth = np.sqrt(abs(cov[1][1]))
    X = np.arange(mu[0]-3*xwidth, mu[0]+3*xwidth,
            xwidth/5.0)
    Y = np.arange(mu[1]-3*ywidth, mu[1]+3*ywidth,
            ywidth/5.0)
    X, Y = np.meshgrid(X, Y)

    # creating a universal function from a python function
    # For details, see http://docs.scipy.org/doc/numpy/reference/ufuncs.html
    ucalc_z = np.frompyfunc(functools.partial(calc_z, mu=mu, cov=cov), 2, 1)

    # calulating z values from corresponding x and y values
    Z = ucalc_z(X, Y)

    ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet, linewidth=1, antialiased=False)
    plt.show()

def main():
    """parses args and passes them to gaussian3d()"""
    if len(sys.argv) < 7:
        print "usage %s <mu0> <mu1> <cov00> <cov01> <cov10> <cov11>"
        return

    # converting argumetns to floating-point values
    argv = map(float, sys.argv[1:])

    # creating mu(average) and covariance
    mu = np.array([[argv[0]], [argv[1]]])
    cov = np.array([[argv[2], argv[3]],[argv[4], argv[5]]])

    gaussian3d(mu, cov)

if __name__ == '__main__':
    main()

以上のコードに適当な名前、例えばgaussian3d.py、を付けて保存する。

$ python gaussian3d.py 0 0 1 0.5 0.5 1

後は上のようなコマンドを実行すれば二次元のガウス分布のグラフが現れる。実際には分布をマウスでグリグリ動かせるのだが、適当にキャプチャした画像を貼っておく。

パラメータを変えながら色々な二次元のガウス分布を描いてみれば、二次元空間までは「納得感」が得られるだろうと期待。