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

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

概要

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

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

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

https://docs.scipy.org/doc/numpy/reference/routines.random.html

　なんとなく「データの生成は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軸グラフの描き方についてはこちらを参考にしました。

[python]matplotlibで左右に２つの軸があるグラフを書く方法 - Qiita

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

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

まとめ

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

静かなる名辞

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

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

概要

方法

実験

まとめ