乱数データと確率密度関数を一緒にplotしてみたかったので、メモ。
概要
scipy.statsでは様々な統計用のユーティリティが提供されています。大抵の分布はあるし、パラメータも好きに設定できます。
numpyにも充実したrandomモジュールがありますが、こちらは分布に従うデータの生成や、データのサンプリングなどしかできません（と思います）。
なんとなく「データの生成はnumpyでできるけど、確率密度関数だとscipy使わないと駄目なのかな？」と思いがちですが、「データの生成」も実はscipyでできるので、numpyを使う必要性はありません。やったね。
方法
この記事では正規分布でやってみます。
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軸グラフの描き方についてはこちらを参考にしました。
結果は、以下のようになりました。
ちゃんとできていますね。
まとめ
scipyでできることがわかった。scipyを信じていれば大抵のことはできる。