In [22]:

    

# カイ二乗分布に従う乱数を指数乱数から生成
import numpy as np
from scipy.stats import uniform, chi2
import matplotlib.pyplot as plt

nbins = 50

# カイ二乗分布のパラメータ（自由度: degrees of freedom）
# この方法では偶数のみ
df = 2
if df % 2 != 0:
    print u"ERROR: 自由度は偶数のみ"
    exit(0)

nu = int(df / 2)

# 逆変換法で一様乱数を指数乱数に変換
# 1つのカイ二乗乱数を得るためにはnu個の指数乱数が必要なので
# 行列化して和を取りやすくする
np.random.seed()
N = 100000
rv = uniform(loc=0.0, scale=1.0)
U = rv.rvs(nu * N).reshape((nu, -1))
X = - np.log(1 - U)

# 指数乱数からカイ二乗乱数を得る
# 各列を足しあわせて1つのカイ二乗乱数を得る
Y = 2 * np.sum(X, axis=0)

# 変換したカイ二乗分布の乱数とpdfを描画
plt.figure(1)
plt.hist(Y, nbins, normed=True)
rv = chi2(df=df)
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((rv.ppf(0.01), rv.ppf(0.99)))

plt.show()


    

In [31]:

    

# ガンマ分布に従う乱数を指数乱数から生成
import numpy as np
from scipy.stats import uniform, gamma
import matplotlib.pyplot as plt

nbins = 50

# ガンマ分布のパラメータ
k = 2.0
theta = 2.0

# 逆変換法で一様乱数を指数乱数に変換
# 1つのガンマ乱数を得るためにはk個の指数乱数が必要なので
# 行列化して和を取りやすくする
np.random.seed()
N = 100000
rv = uniform(loc=0.0, scale=1.0)
U = rv.rvs(k * N).reshape((k, -1))
X = - np.log(1 - U)

# 指数乱数からガンマ乱数を得る
# 各列を足しあわせて1つのガンマ乱数を得る
Y = theta * np.sum(X, axis=0)

# 変換したガンマ分布の乱数とpdfを描画
plt.figure(1)
plt.hist(Y, nbins, normed=True)
rv = gamma(a=k, scale=theta)
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((rv.ppf(0.01), rv.ppf(0.99)))

plt.show()


    

In [41]:

    

# ベータ分布に従う乱数を指数乱数から生成
import numpy as np
from scipy.stats import uniform, beta
import matplotlib.pyplot as plt

nbins = 50

# ベータ分布のパラメータ
a = 5.0
b = 1.0

# 逆変換法で一様乱数を指数乱数に変換
# 1つのガンマ乱数を得るためにはk個の指数乱数が必要なので
# 行列化して和を取りやすくする
np.random.seed()
N = 100000
rv = uniform(loc=0.0, scale=1.0)
U = rv.rvs((a+b) * N).reshape((a+b, -1))
X = - np.log(1 - U)

# 指数乱数からベータ乱数を得る
Y = np.sum(X[0:a,], axis=0) / np.sum(X, axis=0)

# 変換したベータ分布の乱数とpdfを描画
plt.figure(1)
plt.hist(Y, nbins, normed=True)
rv = beta(a=a, b=b)
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((rv.ppf(0.01), rv.ppf(0.99)))

plt.show()