静かなる名辞

pythonとプログラミングのこと



scipyで確率分布のサンプルと確率密度関数を生成する

 乱数データと確率密度関数を一緒にplotしてみたかったので、メモ。

概要

 scipy.statsでは様々な統計用のユーティリティが提供されています。大抵の分布はあるし、パラメータも好きに設定できます。

Statistical functions (scipy.stats) — SciPy v0.16.1 Reference Guide

 numpyにも充実したrandomモジュールがありますが、こちらは分布に従うデータの生成や、データのサンプリングなどしかできません(と思います)。

Random sampling (numpy.random) — NumPy v1.15 Manual

 なんとなく「データの生成はnumpyでできるけど、確率密度関数だとscipy使わないと駄目なのかな?」と思いがちですが、「データの生成」も実はscipyでできるので、numpyを使う必要性はありません。やったね。

方法

 この記事では正規分布でやってみます。

scipy.stats.norm — SciPy v0.16.1 Reference Guide

 scipy.stats.normというものを使うのですが、これは実はクラスに近い使い方ができます。

>>> from scipy import stats
>>> stats.norm
<scipy.stats._continuous_distns.norm_gen object at 0x7f8d839ede10>
>>> type(stats.norm)
<class 'scipy.stats._continuous_distns.norm_gen'>

 scipy.stats._continuous_distns.norm_gen objectというのはわかりづらいですが、早い話がファクトリであり、callするとオブジェクトを返します。

>>> stats.norm()
<scipy.stats._distn_infrastructure.rv_frozen object at 0x7f8da586ef28>

 素性はいまいちわかりませんが、rv_frozen objectという名前から想像できる通り、事実上「確率分布のインスタンス」のように使えます。

 つまり、分布のパラメータを渡してインスタンス化し、インスタンスを使ってそのパラメータの分布に従うサンプルを生成したり、確率密度関数を計算したりできるのですね。これを知ったときは喜びました。

 使えるメソッドはクラスによって微妙に違いがあるようですが、

  • rvs:Random Variates(確率変量)
  • pdf:Probability density function(確率密度関数)

 の2つが今回使いたい「分布に従うデータ」と「確率密度関数」を返すメソッドです。

 なお、これらのメソッドはstats.normのようなクラス(?)から呼ぶことも、インスタンス化したオブジェクト(?)(rv_frozenのオブジェクトと呼ぶべきか……)から呼ぶこともできますが、同じ名前でも使い方が違うので注意が必要です。

>>> from scipy import stats
>>> help(stats.norm.rvs)
Help on method rvs in module scipy.stats._distn_infrastructure:

rvs(*args, **kwds) method of scipy.stats._continuous_distns.norm_gen instance
    Random variates of given type.
    
    Parameters
    ----------
    arg1, arg2, arg3,... : array_like
        The shape parameter(s) for the distribution (see docstring of the
        instance object for more information).
    loc : array_like, optional
        Location parameter (default=0).
    scale : array_like, optional
        Scale parameter (default=1).
    size : int or tuple of ints, optional
        Defining number of random variates (default is 1).
    random_state : None or int or ``np.random.RandomState`` instance, optional
        If int or RandomState, use it for drawing the random variates.
        If None, rely on ``self.random_state``.
        Default is None.
    
    Returns
    -------
    rvs : ndarray or scalar
        Random variates of given `size`.

>>> help(stats.norm().rvs)
Help on method rvs in module scipy.stats._distn_infrastructure:

rvs(size=None, random_state=None) method of scipy.stats._distn_infrastructure.rv_frozen instance

 やりたいことはわかるし、こういう仕様にしたのも理解できるけど、私のような初心者は戸惑う。scipyのドキュメントの簡潔さも相まって。

実験

 こんなコードを動かしてみました。

import numpy as np
import matplotlib.pyplot as plt

from scipy import stats

def main():
    normd = stats.norm(100, 7)
    x = np.arange(70, 130, 0.1)
    pdf = normd.pdf(x)
    samples = normd.rvs(1000)

    fig, ax1 = plt.subplots()
    ax1.hist(samples, bins=30, color="C0")
    ax2 = ax1.twinx()
    ax2.plot(x, pdf, color="C1")
    plt.savefig("result.png")

if __name__ == "__main__":
    main()

 2軸グラフの描き方についてはこちらを参考にしました。

qiita.com


 結果は、以下のようになりました。

result.png
result.png

 ちゃんとできていますね。

まとめ

 scipyでできることがわかった。scipyを信じていれば大抵のことはできる。