matplotlib.mlab.bivariate_normalの使い方

python numpy matplotlib

ML Advent Calendarの記事を書くための実験に二次元ガウス分布の確率密度関数が必要になってmatplotlib.mlab.bivariate_normalの存在を知って使っていたのだけれど，どうもおかしいということに気づいて，引数の意味を勘違いしていたのでメモ．

API Documentを眺めるとこんな感じになっている．

matplotlib.mlab.bivariate_normal(X, Y, sigmax=1.0, sigmay=1.0, mux=0.0, muy=0.0, sigmaxy=0.0)
Bivariate Gaussian distribution for equal shape X, Y.

See bivariate normal at mathworld.

(http://matplotlib.org/api/mlab_api.html より抜粋)

sigmax，sigmay，sigmaxyと書いてあったので分散共分散行列の値を利用して使ってみたけれど，どうもおかしい．自作多変量ガウス分布の密度値と異なる．念のためRのdmvnorm関数で確認してみたところ，どうやら自作関数が正しいようだ．

結局いろいろ確認したところ，引数の意味を勘違いしていたというオチに気づく．分散共分散行列の対角要素はsigmax^2, sigmay^2ということ．すなわち引数に渡す際にsqrtを取ればよかったのだ．引数にはsigmaxyと書かれていたが，正確にはVxyと記述すべきだろう．

使い方サンプルコード

import numpy
import matplotlib
from matplotlib import pyplot
from matplotlib import mlab

mu    = numpy.array( [-2,-2] )
sigma = numpy.array( [[3,1],[1,3]] )

X, Y = numpy.meshgrid( numpy.linspace( -6, 3, 100 ),
                       numpy.linspace( -6, 3, 100 ) )

Z = mlab.bivariate_normal( X, Y,
                           numpy.sqrt(sigma[ 0 ][ 0 ]),
                           numpy.sqrt(sigma[ 1 ][ 1 ]),
                           mu[ 0 ], mu[ 1 ],
                           sigma[ 0 ][ 1 ] )

pyplot.pcolor(X, Y, Z)
pyplot.show()

しかし，関数としても使いづらいので，muをベクトル，sigmaを行列で受け取れるような以下のようなラッパーを自分で書くのがよいと思う．

def 2dgauss_pdf (X, Y, mu, sigma):
    return mlab.bivariate_normal(X, Y,
                                 numpy.sqrt(sigma[ 0 ][ 0 ]),
                                 numpy.sqrt(sigma[ 1 ][ 1 ]),
                                 mu[0], mu[1],
                                 sigma[ 0 ][ 1 ])

近いうちにscipy.statに多変量ガウス分布の確率密度関数が追加される予定らしいという情報もチラっと見たのだけれど，develバージョンのdocumentにもそれらしい記述がなかったのでしばらくは自作関数を利用するのがよいと思った．

自分はこのURLのコードを少し改造したものを自分用のライブラリとして利用することにした．

2014-01-08追記

以下の情報によると SciPy 0.14.0.dev-16fc0a から追加されたとのこと

from scipy.stats import multivariate_normal
var = multivariate_normal(mean=[0,0], cov=[[1,0],[0,1]])
var.pdf([1,0])