SymPyの使い方17 ～統計モジュール - つれづれなる備忘録

　今回はSymPyの統計計算ツールを利用できるStatsモジュールについて紹介する。

1. シンボリックな乱数

　SymPyの統計モジュールであるStatsを利用することで、通常のシンボリック変数にランダム性を加えるシンボリックな乱数を定義することができる。簡単な例として、サイコロを振ったときの目をシンボリックな乱数を定義する関数Dieをsympy.statsからインポートする。

from sympy.stats import Die
X=Die('X',6)
type(X)
>sympy.stats.rv.RandomSymbol

Die('X',6)は6面のサイコロで、変数をXとする。変数Xは単純なシンボリック変数ではなく、type(X)とするとRandomSymbolという型になっていることがわかる。次にXに対してある条件を満たす確率を計算する関数Pをインポートして、サイコロの目が3よりも大きい確率をP(X>3)として求める。

from sympy.stats import P
P(X>3)
>1/2

4,5,6が出る確率なので1/2となる。numpyのような数値計算モジュールだと、ランダム変数を、例えばベクトル列として多数発生させて、条件を満たす乱数の数から確率を計算するが、ぴったり1/2にはならない。このStatsモジュールの特徴としてはサイコロの4,5,6が出る確率が1/2という理論値を計算する点にある。

上の例は6面のサイコロだったが、Die('X2',5)とすると1～5までの目を持つサイコロによるシンボリックな乱数を生成できる。

X2=Die('X2',5)
P(X2>=2)
>4/5

5面のサイコロで2,3,4,5が出る確率は4/5となる。

2. 期待値、分散

　統計の代表的な指標である期待値や分散もStatsモジュールであれば、乱数分布による理論値を簡単に計算することができる。 6面のサイコロの目の期待値、分散はそれぞれStatsモジュールからE,varianceをインポートする。

from sympy.stats import E,variance
E(X)
>7/2
variance(X)
>35/12

E(X),variance(X)とするだけで6面サイコロの期待値、分散の理論値が得られた。シンボリックな乱数Yを追加して、2つの乱数の期待値、分散も求めることができる。

Y=Die('Y',6)
E(X+Y)
>7
variance(X+Y)
>35/6

さらにシンボリック変数aを入れた場合は、期待値や分散の公式通り(E(aX)=aE(X), V(aX)=a²V(X))に理論値が得られる

a=symbols('a')
E(a*X)
>7*a/2
variance(a*X)
>35*a**2/12

3. 正規分布

　よく使われる正規分布による乱数もシンボリックな乱数として扱える。

from sympy.stats import Normal
Z=Normal('Z',0,1)

関数Normalをインポートし、平均0,分散1, 変数名ZとするにはNormal('Z',0,1)とする。期待値や分散をEやvarianceで求めると設定通り期待値0, 分散1が得られる。

E(Z)
>0
variance(Z)
>1

例えば1より大きい値をとる確率は

simplify(P(Z>1))
>-erf(sqrt(2)/2)/2 + 1/2

P(Z>1)でもよいが、冗長なのでsimplifyで簡単化した。またシンボリック変数aを用いることもできる。

P(Z>a)
>erfc(sqrt(2)*a/2)/2

3. まとめ

今回はsympyで統計計算を扱うstatsモジュールと代表的な使い方について紹介した。