In [144]:

    

# 混合分布によってカイ二乗乱数からスチューデントのt分布に従う乱数を生成
# X|y〜N(0, ν/y) && Y〜Xν^2 => X〜t(ν)
import numpy as np
from scipy.stats import chi2, norm, t
import matplotlib.pyplot as plt

# t分布のパラメータ
nu = 10^6

# サンプリング数
N = 100000

# カイ二乗分布に従う乱数を生成
Y = chi2(df=nu).rvs(size=N)

# 生成したカイ二乗分布の乱数をパラメータに持つ正規分布から乱数を生成
# scaleの次元数 = sizeとなっていないとエラーになる
# scaleの各値について1つの乱数を生成してくれる？
# nu=1のときに微妙に理論値とずれるがscipy.statsの乱数生成でも同じ結果なのでOK？
X = norm(loc=0, scale=np.sqrt(nu/Y)).rvs(size=N)

# 上の書き方の意図
# X = np.array([])
# for y in Y:
#     X = np.append(X, norm(loc=0, scale=np.sqrt(nu/y)).rvs(size=1))

# scipy.statsの乱数生成を使ったとき
# X = t(df=nu).rvs(N)

# t分布の裾野は広いので表示範囲を制限するため
X = X[(X>-5) & (X<5)]

plt.figure(1)
nbins = 100
plt.hist(X, nbins, normed=True)
rv = t(df=nu)
x = np.linspace(rv.ppf(0.01), rv.ppf(0.99), 1000)
y = rv.pdf(x)
plt.plot(x, y, 'r-', lw=2)
plt.xlim((-5, 5))

plt.show()


    

In [213]:

    

# 混合分布によってガンマ乱数から負の二項分布に従う乱数を生成
# Y〜Gamma(n,(1-p)/p) && X|y〜Poisson(y) => X〜Neg(n,p)
import numpy as np
from scipy.stats import gamma, poisson, nbinom
import matplotlib.pyplot as plt

# 負の二項分布のパラメータ
n = 6
p = 0.3

# この値まで確率値を計算
K = 50

# サンプリング数
N = 100000

# ガンマ分布に従う乱数を生成
Y = gamma(a=n, scale=(1-p)/p).rvs(size=N)

# 生成したガンマ分布の乱数をパラメータに持つポアソン分布から乱数を生成
X = poisson(mu=Y).rvs(size=N)

# scipy.statsの機能で負の二項分布のパラメータ生成
# X = nbinom(n=n, p=p).rvs(N)

plt.figure(1)
weights = np.ones_like(X) / float(len(X))
nbins = np.arange(-0.5, K, 1.0)
plt.hist(X, nbins, weights=weights)
rv = nbinom(n=n, p=p)
t = np.arange(K)
plt.plot(t, rv.pmf(t), 'r-', lw=2)
plt.xlim((0, K))
plt.show()