Sympy: Creating a Binomial random variable with Poisson number of trials

Sympy: Creating a Binomial random variable with Poisson number of trials – StackOverflow

Sympyにもsympy.statsというのがあるんだ.知らなかった.順を追ってみていく.

単にコバリアンスが欲しい場合,一番簡単なのはnumpy.
(bias=0:unbiased;不偏共分散,bias=1:biased;標本共分散)

import numpy as np

x = np.random.binomial(10, 0.5, 10000)
y = np.random.poisson(10, 10000)
np.cov(x, y)
#X = np.concatenate([[x], [y]])
#np.cov(X)

array([[ 2.42951299, -0.05484492],
[-0.05484492, 9.90795864]])

確率質量関数(離散確率変数なので)を求めたり,
その他ゴニャゴニャしたい時には,scipy.statsが便利(関連[1]).

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import binom, poisson


r1 = binom.rvs(10, 0.5, size=10000)
r2 = poisson.rvs(10, size=10000)
np.cov(r1, r2)
#X = np.concatenate([[r1], [r2]])
#np.cov(X)

array([[ 2.44600744, 0.09936394],
[ 0.09936394, 9.92949295]])

そして,多分それと同様の事(しかも記号的に)ができるんだろうという事で,
sympy.statsについてみていく.

from sympy import Symbol, Eq, simplify
from sympy.stats import P, E, Covariance, Poisson, Binomial
from sympy import init_printing
init_printing()

z1 = Symbol('z1')
x1 = Poisson('x1', 10)
x2 = Binomial('x2', z1, 0.5)
Covariance(x1, x2)#Covariance(x1, x2).evaluate_integral()

Covariance(x1,x2)

期待値を求め,分散を求め,と自分で計算しなくても,
Covariance(記号計算)やcovarianceが用意されている.

z1が整数じゃないので,evalはできない.

代入すれば良い.

Covariance(x1, x2).subs(z1, 10).evaluate_integral()

0

x1 = Poisson('x1', 10)
x2 = Binomial('x2', 10, 0.5)
covariance(x1, x2)

0

独立なケースをみても,いまいちパッとしないが.

質問のケースでは,二項分布の試行回数nがポアソン分布に従う確率変数な場合を考えたい.

Covariance(記号計算)やcovarianceは「subs(z1, x1)」みたいな事ができない.

で,以下の質問にある様に,

How to create a random variable whose parameters are themselves random variables in SymPy? – StackOverflow

z1 = Symbol("z1")
x1 = Poisson("x1", 10)
x2 = Poisson("x2", z1)
Ex2 = E(E(x2).subs(z1, x1))
Vx2 = E(E((x2-Ex2)**2).subs(z1, x1))
cov = E(E((z1-E(x1))*(x2-Ex2)).subs(z1, x1))
print(f'E(x2) = {Ex2}, var(x2) = {Vx2}, cov(x1, x2) = {cov}')

E(x2) = 10, var(x2) = 20, cov(x1, x2) = 10


y1 = np.random.poisson(10, 10000)
y2 = np.random.poisson(y1, 10000)
np.cov(y1, y2)

array([[ 10.10010272, 10.06958472],

[ 10.06958472, 19.89375994]])

はできるけど,

z1 = Symbol("z1")
x1 = Poisson("x1", 10)
x2 = Binomial("x2", z1, 0.5)
Ex2 = E(E(x2).subs(z1, x1))
Vx2 = E(E((x2-Ex2)**2).subs(z1, x1))
cov = E(E((z1-E(x1))*(x2-Ex2)).subs(z1, x1))
print(f'E(x2) = {Ex2}, var(x2) = {Vx2}, cov(x1, x2) = {cov}')

ValueError: z1 is not an integer

はできない.

Source code for sympy.stats.drv_types

Source code for sympy.stats.frv_types

nは整数じゃないと駄目っぽい.

ちなみに,期待される答えをNumpyやscipy.statsで求めるなら,

y3 = np.random.poisson(10, 100000)
y4 = np.random.binomial(y3, 0.5, 100000)
np.cov(y3, y4)

array([[ 10.01005157, 5.00417544],
[ 5.00417544, 4.99903657]])

r3 = poisson.rvs(10, size=100000)
r4 = binom.rvs(r3, 0.5, size=100000)
np.cov(r3, r4)

array([[ 9.99372993, 5.00735952],
[ 5.00735952, 5.0200852 ]])

 
 
 
関連:
[1] 正規分布の面積~累積確率(累積分布関数と生存関数)の話

広告
カテゴリー: 未分類 パーマリンク

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中