scipy.stat.gaussian_kde()による確率密度推定 - (主に)プログラミングのメモ

scipy.stats のインポート：

>>> from scipy import stats

Gaussian kernel density estimation のためのクラスの定義：

>>> class GaussianKernelDensityEstimation(object):
    """docstring for GaussianKernelDensityEstimation"""
    def __init__(self):
        self.gkde = None

    def fit(self, X):
        """docstring for fit"""
        self.gkde = stats.gaussian_kde(X.T)
        return self

    def predict(self, X):
        return self.gkde.evaluate(X.T)

２次元空間上に乱数を生成：

>>> X = np.random.randn(400,2)    # 標準正規分布 N(0,1) に従う乱数を400個 ＆ ２次元
>>> print X   # サイズ 400 x 2 
[[ -2.58007839e-01   5.48065715e-01]
 [  1.38392260e+00   8.88118290e-01]
 [ -5.54353637e-01  -4.50188476e-01]
 [  1.50815543e+00  -8.12178615e-02]
 [  6.12684130e-02  -8.58764228e-01]
............................
............................
 [ -6.05231082e-01   6.09153545e-01]
 [  8.76175253e-02   1.27423251e+00]
 [ -2.50029715e+00   4.42584247e-01]
 [  7.49176170e-02   1.46820415e+00]]

確率密度の値を得たい点の集合を作る：

>>> xmin, ymin = X.min(axis=0)
>>> xmax, ymax = X.max(axis=0)
>>> xx = np.linspace(xmin - 0.5, xmax + 0.5, 32)
>>> yy = np.linspace(ymin - 0.5, ymax + 0.5, 32)
>>> xg, yg = np.meshgrid(xx, yy)
>>> grid_coords = np.c_[xg.ravel(), yg.ravel()]
>>> print grid_coords
[[-4.00812329 -2.66008256]
 [-3.95031732 -2.66008256]
 [-3.89251135 -2.66008256]
 ..., 
 [ 3.21762312  3.69041801]
 [ 3.27542909  3.69041801]
 [ 3.33323506  3.69041801]]

確率密度関数の推定：

>> kde = GaussianKernelDensityEstimation()   # gaussian kde のオブジェクトを作って
>>> kde.fit(X)    # ここで確率密度を推定している。
<__main__.GaussianKernelDensityEstimation object at 0xa16042c>

grid_coords に含まれる個々の点における確率密度を求める：

>>> zz = kde.predict(grid_coords)

>>> len(grid_coords)   # ２次元空間上の点が 1024
1024
>>> len(zz)  #当然，それに対応する確率密度の個数も 1024
1024
>>> 
>>> for i in range(len(grid_coords)):
    print "%f,%f,%f" %(grid_coords[i][0], grid_coords[i][1], zz[i])

    
-4.008123,-2.660083,0.000002
-3.771305,-2.660083,0.000004
-3.534487,-2.660083,0.000006
....................
....................
2.859599,3.690418,0.000000
3.096417,3.690418,0.000000
3.333235,3.690418,0.000000
>>>

上記の結果を３列目(=確率密度)でソートした結果：

0.254601,-0.201824,0.149000
0.254601,0.003031,0.146694
0.017783,-0.201824,0.146501
0.017783,0.003031,0.144133
....................
....................
-3.060851,3.280708,0.000000
-3.060851,3.075853,0.000000
-2.824033,3.690418,0.000000

原点周辺での確率密度が大きくなっているようだ。

原点における確率密度を求めると

>>> Y = np.array([[0.0,0.0]])
>>> Y
array([[ 0.,  0.]])
>>> kde.predict(Y)
array([ 0.14375202])

サンプル数(= X のサイズ）が少ないため，原点に近いからといって最大となっていない。

ちょっと悔しいから，
50,000個のサンプルを与えて原点およびその周辺の点における確率密度を求めた。
原点から遠ざかるほど確率密度が小さくなり，原点から等距離にある４つの点での密度がほぼ等しく収束していることが分かる。

>>> X = np.random.randn(50000,2)
>>> kde = GaussianKernelDensityEstimation()
>>> kde.fit(X)
<__main__.GaussianKernelDensityEstimation object at 0xa16b22c>
>>> Y = np.array([[0.0,0.0]])
>>> kde.predict(Y)
array([ 0.14903213])
>>> Y = np.array([[1,0], [-1,0], [0,1], [0,-1]])
>>> kde.predict(Y)
array([ 0.09236879,  0.09647149,  0.09689293,  0.09497972])
>>> Y = np.array([[2,0], [-2,0], [0,2], [0,-2]])
>>> kde.predict(Y)
array([ 0.0223722 ,  0.02150799,  0.0207469 ,  0.02212046])