SymPyの使い方19 ～統計モジュール3

　前回のシンボリックな乱数から乱数生成、任意の確率分布に従うシンボリックな乱数に関して、今回は度数分布を表示する方法について紹介する。

atatat.hatenablog.com

1. 乱数の生成(変更)
2. 生成された乱数の度数分布表示
3. まとめ

1. 乱数の生成(変更)

　前回紹介したシンボリックな乱数から乱数を生成する方法に関して、恐らくgoogle colabのpythonバージョンのアップデートがあったため変更が発生している。前回と同じく６面サイコロのシンボリックな乱数を生成するには以下のようにする　

from sympy import *
from sympy.stats import Die, sample
X=Die('X',6)

前回は単にsample(X)とすれば乱数が発生していたが、sample(X)がiterator objectに扱いが変更になったため、以下のようにして乱数を生成する。

next(sample(X))
>1

なおsample関数を使用するたびに以下の警告メッセージが表記されるが、内容としてはsampleの返り値がversion 1.7からiterator objectに変更されたことが書かれている。

/usr/local/lib/python3.7/dist-packages/sympy/stats/rv.py:1104: UserWarning: 
The return type of sample has been changed to return an iterator
object since version 1.7. For more information see
https://github.com/sympy/sympy/issues/19061
  warnings.warn(filldedent(message))

2. 生成された乱数の度数分布表示

　シンボリックな乱数から乱数を生成し、ヒストグラムを用いて度数分布を表示する。6面サイコロの乱数Xは以下のようにforループで繰り返し乱数を1000回発生させて、リストxxに追加する。ヒストグラムを生成するにはhist()を用いるが、リストxxをnumpyアレイに変換しておく。

import matplotlib.pyplot as plt
import numpy as np
xx=[]
for i in np.arange(1000):
 xx.append(next(sample(X)))
plt.hist(np.array(xx),bins=11)

以下のように1から6までほぼ同じ頻度(若干6が少ない)で乱数が発生していることがわかる。

6面サイコロと同様の手順で平均0,分散1の正規分布のシンボリックな乱数から度数分布を表示する。

from sympy.stats import Normal
Z=Normal('Z',0,1)
zz=[]
for i in np.arange(1000):
  zz.append(next(sample(Z)))
plt.hist(np.array(zz),bins=20)

平均0を中心としたガウシアン形状の頻度分布になっていることが確認できる。

任意の確率分布に従ったシンボリックな乱数の度数分布も表示することができるが、あらかじめ用意されているDieやNormalよりもfor文の処理が遅いことに注意が必要。

from sympy import exp, Symbol, Interval, oo
from sympy.stats import ContinuousRV
x=Symbol('x')
pdf=exp(-x)
Z2 = ContinuousRV(x, pdf, set=Interval(0, oo))
zz2=[]
for i in np.arange(1000):
 zz2.append(next(sample(Z2)))
plt.hist(np.array(zz2),bins=20)

任意の離散的な分布についても同様に度数分布を表示できる。for文の処理が遅い点は連続分布と同様。

from sympy import S
from sympy.stats import DiscreteRV
p=S(1)/2
x = Symbol('x', integer=True, positive=True)
pdf = p*(1 - p)**(x - 1)
D = DiscreteRV(x, pdf, set=S.Naturals)
d=[]
for i in np.arange(1000):
 d.append(next(sample(D)))
plt.hist(np.array(d),bins=20)