ホーム » 分布
「分布」カテゴリーアーカイブ
二項分布と超幾何分布
二項分布
- 抽出方法: 復元抽出
- 試行の独立性: 各試行は独立
- 確率: 各試行で成功確率は一定
超幾何分布
- 抽出方法: 非復元抽出
- 試行の独立性: 各試行は独立していない
- 確率: 各試行で成功確率が変化
グラフの比較
今回の実験では、二項分布と超幾何分布の違いを視覚的に比較します。二項分布は、復元抽出を行う場合に適用され、各試行が独立しており、成功確率が一定です。一方、超幾何分布は、非復元抽出を行う場合に適用され、試行ごとに成功確率が変化します。実験では、母集団のサイズ、成功対象の数、抽出回数を変え、それぞれのグラフを描画します。特に、サンプルサイズが大きい場合や母集団が小さい場合など、分布の形状に顕著な違いが現れる条件に注目して比較を行います。これにより、抽出方法による分布の違いを視覚的に確認し、二項分布と超幾何分布の特性を理解します。
復元抽出と非復元抽出の比較(広範囲のパラメータ設定)
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binom, hypergeom
# グラフを描画する関数
def plot_comparison(M, N_A, n, ax):
# 非復元抽出(超幾何分布)
rv_hypergeom = hypergeom(M, N_A, n)
x_values = np.arange(0, n+1)
pmf_hypergeom = rv_hypergeom.pmf(x_values)
# 復元抽出(二項分布)
p = N_A / M # 豆Aを引く確率
rv_binom = binom(n, p)
pmf_binom = rv_binom.pmf(x_values)
# グラフ描画
ax.plot(x_values, pmf_hypergeom, 'bo-', label='非復元抽出 (超幾何分布)', markersize=5)
ax.plot(x_values, pmf_binom, 'ro-', label='復元抽出 (二項分布)', markersize=5)
condition_text = f'M = {M}, N_A = {N_A}, n = {n}'
ax.text(0.05, 0.95, condition_text, transform=ax.transAxes,
fontsize=12, verticalalignment='top')
ax.set_xlabel('豆Aの数')
ax.set_ylabel('確率')
ax.grid(True)
ax.legend()
# パラメータ比を一定に保ったグラフを作成
M_values = [50, 150, 250, 350, 450]
N_A_values = [15, 45, 75, 105, 135] # M の 30%
n_values = [7, 22, 37, 52, 67] # M の 15%
# サブプロットを作成(縦に並べる)
fig, axs = plt.subplots(5, 1, figsize=(10, 30))
# 異なるパラメータ設定でグラフを描画
for i in range(5):
plot_comparison(M=M_values[i], N_A=N_A_values[i], n=n_values[i], ax=axs[i])
# レイアウトの自動調整
plt.tight_layout()
# 余白の手動調整
plt.subplots_adjust(right=0.95)
plt.show()
復元抽出と非復元抽出の比較(顕著な違いが現れる条件)
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binom, hypergeom
# グラフを描画する関数
def plot_comparison(M, N_A, n, ax):
# 非復元抽出(超幾何分布)
rv_hypergeom = hypergeom(M, N_A, n)
x_values = np.arange(0, n+1)
pmf_hypergeom = rv_hypergeom.pmf(x_values)
# 復元抽出(二項分布)
p = N_A / M # 豆Aを引く確率
rv_binom = binom(n, p)
pmf_binom = rv_binom.pmf(x_values)
# グラフ描画
ax.plot(x_values, pmf_hypergeom, 'bo-', label='非復元抽出 (超幾何分布)', markersize=5)
ax.plot(x_values, pmf_binom, 'ro-', label='復元抽出 (二項分布)', markersize=5)
condition_text = f'M = {M}, N_A = {N_A}, n = {n}'
ax.text(0.05, 0.95, condition_text, transform=ax.transAxes,
fontsize=12, verticalalignment='top')
ax.set_xlabel('豆Aの数')
ax.set_ylabel('確率')
ax.grid(True)
ax.legend()
# サブプロットを作成(縦に並べる)
fig, axs = plt.subplots(3, 1, figsize=(10, 18)) # 高さを調整
# 条件1: サンプルサイズが大きい場合
plot_comparison(M=100, N_A=30, n=80, ax=axs[0])
# 条件2: 全体のサイズが小さい場合
plot_comparison(M=20, N_A=6, n=15, ax=axs[1])
# 条件3: 極端な割合の場合
plot_comparison(M=100, N_A=90, n=30, ax=axs[2])
# レイアウトの自動調整
plt.tight_layout()
# 余白の手動調整
plt.subplots_adjust(right=0.95)
plt.show()
考察
二項分布と超幾何分布の違いが顕著になるのは、サンプルサイズが全体に対して大きい場合や、母集団が小さい場合です。例えば、全体の豆の数が100個で80個を抽出する場合や、母集団が20個でそのうち15個を抽出する場合、非復元抽出では次の試行に与える影響が大きくなり、復元抽出との差が明確に現れます。また、成功確率が極端に高いか低い場合も、非復元抽出の影響で分布が変わりやすくなります。一方、サンプルサイズが小さい場合や母集団が非常に大きい場合には、各抽出の影響が少なく、二項分布と超幾何分布の違いは目立たなくなります。
混合ガウスモデルによる観測データの分布特性推定
概要:
本実験では、観測データが複数の異なる分布の混合によって生成されている可能性を考慮し、その分布特性を推定する手法を検証する。具体的には、2つの正規分布を混合したデータを生成し、EMアルゴリズムを用いてその混合ガウスモデルをフィッティングする。実験の結果、生成されたデータの背後にある各分布の平均、標準偏差、および混合割合を高精度に推定できることを確認した。本手法は、複雑な観測データに対する分布特性の解析やデータの要因分解に有効である。
2つの成分
import numpy as np
from sklearn.mixture import GaussianMixture
import matplotlib.pyplot as plt
# 1. データ生成: 混合正規分布から乱数を生成
np.random.seed(0)
n_samples = 1000
# 2つの正規分布のパラメータ設定
mean1, std1 = 0, 1 # 正規分布1
mean2, std2 = 5, 1.5 # 正規分布2
weight1, weight2 = 0.4, 0.6 # 混合割合
# 2つの分布からデータ生成
data1 = np.random.normal(mean1, std1, int(weight1 * n_samples))
data2 = np.random.normal(mean2, std2, int(weight2 * n_samples))
# 混合データの作成
data = np.hstack([data1, data2])
np.random.shuffle(data)
# 2. 混合モデルの構築とフィッティング (EMアルゴリズムを使用)
gmm = GaussianMixture(n_components=2, random_state=0)
gmm.fit(data.reshape(-1, 1))
# 推定されたパラメータ
means = gmm.means_.flatten()
std_devs = np.sqrt(gmm.covariances_).flatten()
weights = gmm.weights_
print("推定された平均:", means)
print("推定された標準偏差:", std_devs)
print("推定された混合割合:", weights)
# 3. 結果のプロット
x = np.linspace(-3, 10, 1000)
pdf1 = weights[0] * (1 / (std_devs[0] * np.sqrt(2 * np.pi)) * np.exp(-(x - means[0])**2 / (2 * std_devs[0]**2)))
pdf2 = weights[1] * (1 / (std_devs[1] * np.sqrt(2 * np.pi)) * np.exp(-(x - means[1])**2 / (2 * std_devs[1]**2)))
plt.hist(data, bins=30, density=True, alpha=0.5, color='gray', label='データのヒストグラム')
plt.plot(x, pdf1 + pdf2, 'r-', lw=2, label='フィッティングされたモデル')
plt.plot(x, pdf1, 'b--', lw=2, label='成分 1')
plt.plot(x, pdf2, 'g--', lw=2, label='成分 2')
plt.title('データとフィッティングされた混合ガウスモデル')
plt.xlabel('値')
plt.ylabel('密度')
plt.legend()
plt.show()
3つの成分
import numpy as np
from sklearn.mixture import GaussianMixture
import matplotlib.pyplot as plt
# 1. データ生成: 混合正規分布から乱数を生成
np.random.seed(0)
n_samples = 1000
# 3つの正規分布のパラメータ設定
mean1, std1 = 0, 1 # 正規分布1
mean2, std2 = 5, 1.5 # 正規分布2
mean3, std3 = 10, 2 # 正規分布3
weight1, weight2, weight3 = 0.3, 0.4, 0.3 # 混合割合
# 3つの分布からデータ生成
data1 = np.random.normal(mean1, std1, int(weight1 * n_samples))
data2 = np.random.normal(mean2, std2, int(weight2 * n_samples))
data3 = np.random.normal(mean3, std3, int(weight3 * n_samples))
# 混合データの作成
data = np.hstack([data1, data2, data3])
np.random.shuffle(data)
# 2. 混合モデルの構築とフィッティング (EMアルゴリズムを使用)
gmm = GaussianMixture(n_components=3, random_state=0)
gmm.fit(data.reshape(-1, 1))
# 推定されたパラメータ
means = gmm.means_.flatten()
std_devs = np.sqrt(gmm.covariances_).flatten()
weights = gmm.weights_
print("推定された平均:", means)
print("推定された標準偏差:", std_devs)
print("推定された混合割合:", weights)
# 3. 結果のプロット
x = np.linspace(-5, 15, 1000)
pdf1 = weights[0] * (1 / (std_devs[0] * np.sqrt(2 * np.pi)) * np.exp(-(x - means[0])**2 / (2 * std_devs[0]**2)))
pdf2 = weights[1] * (1 / (std_devs[1] * np.sqrt(2 * np.pi)) * np.exp(-(x - means[1])**2 / (2 * std_devs[1]**2)))
pdf3 = weights[2] * (1 / (std_devs[2] * np.sqrt(2 * np.pi)) * np.exp(-(x - means[2])**2 / (2 * std_devs[2]**2)))
plt.hist(data, bins=30, density=True, alpha=0.5, color='gray', label='データのヒストグラム')
plt.plot(x, pdf1 + pdf2 + pdf3, 'r-', lw=2, label='フィッティングされたモデル')
plt.plot(x, pdf1, 'b--', lw=2, label='成分 1')
plt.plot(x, pdf2, 'g--', lw=2, label='成分 2')
plt.plot(x, pdf3, 'y--', lw=2, label='成分 3')
plt.title('データとフィッティングされた混合ガウスモデル')
plt.xlabel('値')
plt.ylabel('密度')
plt.legend()
plt.show()
2019 Q3(1)
一様分布に於いて最大値が十分統計量であることを示しました。
コード
一様分布の最大値Yを元に、不偏推定量を計算します。
import numpy as np
import matplotlib.pyplot as plt
# シミュレーションパラメータ
n = 10 # サンプルサイズ
theta_true = 10 # 真のθ
num_simulations = 1000 # シミュレーション回数
# θの推定値を格納するリスト
theta_estimates = []
# シミュレーション開始
for _ in range(num_simulations):
# (0, theta_true) の一様分布からn個のサンプルを生成
samples = np.random.uniform(0, theta_true, n)
# 最大値 Y = max(X_1, ..., X_n)
Y = np.max(samples)
# Y から θ を推定 (最大値を使って推定)
theta_hat = Y * (n / (n - 1)) # 推定式: Y * n / (n-1)
# 推定されたθをリストに保存
theta_estimates.append(theta_hat)
# θの推定値の平均を計算
theta_mean = np.mean(theta_estimates)
# 推定結果をヒストグラムで表示
plt.hist(theta_estimates, bins=30, alpha=0.7, color='blue', edgecolor='black', label='推定されたθ')
plt.axvline(x=theta_true, color='red', linestyle='--', label=f'真のθ = {theta_true}')
plt.axvline(x=theta_mean, color='green', linestyle='--', label=f'推定されたθの平均 = {theta_mean:.2f}')
plt.title(f'Y (最大値) に基づく θ の推定値のヒストグラム ({n} サンプル)')
plt.xlabel('推定された θ')
plt.ylabel('頻度')
plt.legend()
plt.show()
# 推定結果をテキストで表示
print(f'真のθ: {theta_true}')
print(f'推定されたθの平均: {theta_mean:.2f}')
真のθ: 10
推定されたθの平均: 10.12
推定されたθは真のθに近い値を取りました。最大値Yだけを用いてθが正確に推定できたので、Yはθに関する十分統計量であることが確認できました。
2019 Q2(4)
損失関数の期待値を導出しそれが最小となるパラメータαを求める問題をやりました。
コード
数値シミュレーションによりR(α)が最小値になるαを見つけます。
# 2019 Q2(4) 2024.9.16
import numpy as np
import matplotlib.pyplot as plt
# 損失関数 R(α) の定義
def R(alpha):
return alpha + 2 / alpha - 2
# α の範囲を設定 (0.1 から 3 まで、細かいステップで計算)
alpha_range = np.linspace(0.1, 3, 1000)
# 損失関数 R(α) の計算
R_values = R(alpha_range)
# 最小値を取る α の確認 (numpy の argmin を使って最小値を見つける)
min_index = np.argmin(R_values)
min_alpha = alpha_range[min_index]
min_R = R_values[min_index]
# グラフのプロット
plt.figure(figsize=(8, 6))
plt.plot(alpha_range, R_values, label='R(α)')
plt.axvline(np.sqrt(2), color='red', linestyle='--', label=r'$\alpha = \sqrt{2}$')
plt.scatter(min_alpha, min_R, color='red', zorder=5, label=f'最小値: α={min_alpha:.2f}, R(α)={min_R:.2f}')
plt.xlabel(r'$\alpha$')
plt.ylabel(r'$R(\alpha)$')
plt.title(r'$R(\alpha)$ の最小値確認 (数値シミュレーション)')
plt.legend()
plt.grid(True)
plt.show()
# 最小値の表示
print(f'数値シミュレーションでの最小値 α: {min_alpha}, R(α): {min_R}')
数値シミュレーションでの最小値 α: 1.4150150150150151, R(α): 0.8284275786822475
手計算と近い値になりました。
極値付近が少し平らなので、αの範囲を1.0~2.0に変更し拡大してみます。
import numpy as np
import matplotlib.pyplot as plt
# 損失関数 R(α) の定義
def R(alpha):
return alpha + 2 / alpha - 2
# α の範囲を 1 から 2 に設定
alpha_range = np.linspace(1, 2, 1000)
# 損失関数 R(α) の計算
R_values = R(alpha_range)
# 最小値を取る α の確認 (numpy の argmin を使って最小値を見つける)
min_index = np.argmin(R_values)
min_alpha = alpha_range[min_index]
min_R = R_values[min_index]
# グラフのプロット
plt.figure(figsize=(8, 6))
plt.plot(alpha_range, R_values, label='R(α)')
plt.axvline(np.sqrt(2), color='red', linestyle='--', label=r'$\alpha = \sqrt{2}$')
plt.scatter(min_alpha, min_R, color='red', zorder=5, label=f'最小値: α={min_alpha:.4f}, R(α)={min_R:.4f}')
plt.xlabel(r'$\alpha$')
plt.ylabel(r'$R(\alpha)$')
plt.title(r'$R(\alpha)$ の最小値確認 (αの範囲: 1 から 2)')
plt.legend()
plt.grid(True)
plt.show()
# 最小値の表示
print(f'数値シミュレーションでの最小値 α: {min_alpha}, R(α): {min_R}')
数値シミュレーションでの最小値 α: 1.4144144144144144, R(α): 0.8284271532679175
最小値が取れていることが分かりました。
2019 Q2(3)
(2)で求めた確率密度関数において確率変数Uの逆数の期待値を求めました。
コード
数式を使って1/Uの期待値を求めます。
# 2019 Q2(3) 2024.9.15
import sympy as sp
# 変数の定義
u, lambda_value = sp.symbols('u lambda', positive=True, real=True)
# ガンマ分布の確率密度関数 g(u) = λ^2 * u * exp(-λ * u)
g_u = lambda_value**2 * u * sp.exp(-lambda_value * u)
# 1/U の期待値の積分
integrand = (1/u) * g_u # 1/u * g(u)
# 積分の実行
expectation = sp.integrate(integrand, (u, 0, sp.oo))
# 結果の簡略化
expectation_simplified = sp.simplify(expectation)
# 結果の表示
display(expectation_simplified)
手計算と一致します。
次に、数値シミュレーションを行います。
import numpy as np
import matplotlib.pyplot as plt
# パラメータ設定
lambda_value = 2 # λ=2
num_samples = 100000 # サンプルの数
# 指数分布に従う乱数 X1, X2 を生成
X1 = np.random.exponential(scale=1/lambda_value, size=num_samples)
X2 = np.random.exponential(scale=1/lambda_value, size=num_samples)
# U = X1 + X2 の計算
U = X1 + X2
# 1/U の計算
inv_U = 1 / U
# シミュレーションの結果からの平均値
simulation_result = np.mean(inv_U)
# 理論値 E[1/U] = λ
theoretical_value = lambda_value
# 結果の表示
print(f'シミュレーションによる E[1/U]: {simulation_result}')
print(f'理論値 E[1/U]: {theoretical_value}')
# ヒストグラムの描画
plt.hist(inv_U, bins=50, density=True, alpha=0.7, label='シミュレーション結果')
plt.axvline(theoretical_value, color='r', linestyle='--', label=f'理論値 E[1/U] = {theoretical_value}')
plt.xlabel('1/U')
plt.ylabel('確率密度')
plt.title('シミュレーションによる 1/U の分布と理論値の比較')
plt.legend()
plt.show()
シミュレーションによる E[1/U]: 1.9928303711222939
理論値 E[1/U]: 2
シミュレーションによりE[1/U]は理論値に近づきました。分散は発散するためヒストグラムにすると大きなバラつきが見られます。
2019 Q2(2)
独立な指数分布に従う2変数の和の確率密度関数を求めました。
コード
数式を使ってUの確率密度関数を求めます。
# 2019 Q2(2) 2024.9.14
import sympy as sp
# 変数の定義
x, u, lambda_value = sp.symbols('x u lambda', positive=True, real=True)
# 指数分布の確率密度関数 (PDF) f(x)
f_x = lambda_value * sp.exp(-lambda_value * x)
# 畳み込み積分を計算して、g(u) = f(x) * f(u - x) の形式にする
g_u = sp.integrate(f_x * lambda_value * sp.exp(-lambda_value * (u - x)), (x, 0, u))
# 簡略化
g_u_simplified = sp.simplify(g_u)
# 結果の表示
display(g_u_simplified)
X1 (X2)と、Uの確率密度関数を重ねて描画してみます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma
# パラメータ設定
lambda_value = 2 # λ=2
x_range = np.linspace(0, 5, 1000) # X1, X2 の範囲
# X1, X2 の指数分布の確率密度関数 (PDF) - 同じなので1つにまとめる
pdf_X1_X2 = lambda_value * np.exp(-lambda_value * x_range)
# U = X1 + X2 は形状パラメータk=2、スケールパラメータθ=1/λのガンマ分布
u_range = np.linspace(0, 10, 1000)
pdf_U = gamma.pdf(u_range, a=2, scale=1/lambda_value)
# グラフのプロット
plt.figure(figsize=(8, 6))
# X1 (X2) の指数分布の描画
plt.plot(x_range, pdf_X1_X2, label='X1 (X2) (指数分布)', color='blue', linestyle='--')
# U のガンマ分布の描画
plt.plot(u_range, pdf_U, label='U = X1 + X2 (ガンマ分布)', color='red')
# グラフの装飾
plt.title('X1 (X2) と U の確率密度関数の比較')
plt.xlabel('値')
plt.ylabel('確率密度')
plt.legend()
# グラフの表示
plt.show()
グラフの形状から、U = X1 + X2の関係は想像しづらいですね。
ガンマ分布の形状パラメータを1~2に変化させて、X1 (X2)からUへの変化を確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma
# パラメータ設定
lambda_value = 2 # λ=2
x_range = np.linspace(0, 5, 1000) # X1, X2, U の範囲を0~5に設定
# ガンマ分布の形状パラメータを1から2まで0.1ずつ変化させる
shape_params = np.arange(1, 2.1, 0.1) # 1から2までの形状パラメータを0.1ステップで変化
# グラフのプロット
plt.figure(figsize=(8, 6))
# 形状パラメータが1から2に変化するガンマ分布の描画
for k in shape_params:
gamma_pdf = gamma.pdf(x_range, a=k, scale=1/lambda_value)
plt.plot(x_range, gamma_pdf, label=f'形状パラメータ k={k:.1f}')
# X1 (X2) の指数分布の描画 - 同じなので1つにまとめる
pdf_X1_X2 = lambda_value * np.exp(-lambda_value * x_range)
plt.plot(x_range, pdf_X1_X2, label='X1 (X2) (指数分布)', color='blue', linestyle='--')
# グラフの装飾
plt.title('指数分布からガンマ分布への変化 (形状パラメータの0.1ステップ変化)')
plt.xlabel('値')
plt.ylabel('確率密度')
plt.legend()
# グラフの表示
plt.show()
X1 (X2)からUへの変化が可視化できました。
次に、X1 (X2)と、Uの累積分布関数を重ねて描画してみます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import expon, gamma
# パラメータ設定
lambda_value = 2 # λ=2
x_range = np.linspace(0, 5, 1000) # X1, X2 の範囲
# X1, X2 の指数分布の累積分布関数 (CDF) - 同じなので1つにまとめる
cdf_X1_X2 = expon.cdf(x_range, scale=1/lambda_value)
# U = X1 + X2 は形状パラメータk=2、スケールパラメータθ=1/λのガンマ分布
u_range = np.linspace(0, 10, 1000)
cdf_U = gamma.cdf(u_range, a=2, scale=1/lambda_value)
# グラフのプロット
plt.figure(figsize=(8, 6))
# X1 (X2) の指数分布のCDFの描画
plt.plot(x_range, cdf_X1_X2, label='X1 (X2) (指数分布 CDF)', color='blue', linestyle='--')
# U のガンマ分布のCDFの描画
plt.plot(u_range, cdf_U, label='U = X1 + X2 (ガンマ分布 CDF)', color='red')
# グラフの装飾
plt.title('X1 (X2) と U の累積分布関数 (CDF) の比較')
plt.xlabel('値')
plt.ylabel('累積確率')
plt.legend()
# グラフの表示
plt.show()
グラフの形状から、U = X1 + X2の関係を想像しやすくなりました。
2019 Q2(1)
独立な指数分布に従う2変数の和の期待値を求めました。
コード
数式を使ってUの期待値を求めます
import sympy as sp
# 変数の定義
x = sp.Symbol('x', positive=True)
lambda_value = sp.Symbol('lambda', positive=True)
# 指数分布の確率密度関数 (PDF)
f_x = lambda_value * sp.exp(-lambda_value * x)
# 期待値の式を計算 (積分)
E_X1 = sp.integrate(x * f_x, (x, 0, sp.oo))
# 期待値
E_X1_simplified = sp.simplify(E_X1)
# X1 + X2 の場合 (独立性より)
E_U = 2 * E_X1_simplified
# 結果の表示
display(E_U)
手計算と一致します
次にUに従う乱数を生成し分布を確認します。また形状パラメータ2のガンマ分布と重ねて描画します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma
# パラメータ設定
lambda_value = 2 # 例としてλ=2
num_samples = 1000 # サンプルの数
# 指数分布に従う乱数の生成
X1 = np.random.exponential(scale=1/lambda_value, size=num_samples)
X2 = np.random.exponential(scale=1/lambda_value, size=num_samples)
# U = X1 + X2 の計算
U = X1 + X2
# 理論的な期待値
expected_value = 2 / lambda_value
# シミュレーション結果のヒストグラムのプロット (確率密度にスケール)
plt.hist(U, bins=30, density=True, alpha=0.7, label='シミュレーション結果')
# ガンマ分布 (形状パラメータ k=2, スケールパラメータ θ=1/λ) の確率密度関数 (PDF) を計算
k = 2 # 形状パラメータ (k = 2)
theta = 1 / lambda_value # スケールパラメータ (θ = 1/λ)
# ガンマ分布に基づくx軸の範囲
x = np.linspace(0, max(U), 1000)
# ガンマ分布の確率密度関数 (PDF) を計算
gamma_pdf = gamma.pdf(x, a=k, scale=theta)
# ガンマ分布の折れ線グラフをプロット
plt.plot(x, gamma_pdf, 'r-', label=f'ガンマ分布 (k={k}, θ={theta:.2f})')
# グラフの装飾
plt.axvline(expected_value, color='g', linestyle='--', label=f'理論値 = {expected_value:.2f}')
plt.xlabel('U = X1 + X2')
plt.ylabel('確率密度')
plt.title('Uの分布とガンマ分布の比較')
plt.legend()
plt.show()
# 実際のサンプル平均値と理論的期待値の比較
print(f'理論的な期待値: {expected_value}')
print(f'シミュレーションからの平均値: {np.mean(U)}')
理論的な期待値: 1.0
シミュレーションからの平均値: 0.9902440726761544
Uの分布と形状パラメータ2のガンマ分布は重なりました。Uはガンマ分布に従うことが分かります。
2021 Q3(1)
ポアソン分布のモーメント母関数を求めました。
コード
数式を使った計算
# 2021 Q3(1) 2024.8.25
import sympy as sp
# 変数の定義
s, x, lambda_ = sp.symbols('s x lambda_', real=True, positive=True)
# ポアソン分布の確率質量関数 (PMF)
poisson_pmf = (lambda_**x * sp.exp(-lambda_)) / sp.factorial(x)
# モーメント母関数 M_X(s) の定義
M_X = sp.summation(sp.exp(s*x) * poisson_pmf, (x, 0, sp.oo))
# 結果を簡略化
M_X_simplified = sp.simplify(M_X)
# 結果を表示
M_X_simplified
無限級数がe^xの形に変換されないようです。
exp(s) を t に置き換えて、最後に戻してみます。
# 2021 Q3(1) 2024.8.25
import sympy as sp
# 変数の定義
s, x, lambda_, t = sp.symbols('s x lambda_ t', real=True, positive=True)
# exp(s) を t に置き換え
summand = (lambda_ * t)**x / sp.factorial(x)
# モーメント母関数 M_X(s) の定義
M_X = sp.exp(-lambda_) * sp.summation(summand, (x, 0, sp.oo))
# 簡略化
M_X_simplified = sp.simplify(M_X)
# 最後に t を exp(s) に戻す
M_X_final = M_X_simplified.subs(t, sp.exp(s))
# 結果を表示
M_X_final
手計算と同じ形になりました。
モーメント母関数を使って期待値と分散を求めます。
import sympy as sp
# 変数の定義
s, lambda_ = sp.symbols('s lambda_', real=True, positive=True)
# モーメント母関数 M_X(s) の定義
M_X = sp.exp(lambda_ * (sp.exp(s) - 1))
# 1次モーメント(期待値)の計算
M_X_prime = sp.diff(M_X, s)
expectation = M_X_prime.subs(s, 0)
# 2次モーメントの計算
M_X_double_prime = sp.diff(M_X_prime, s)
second_moment = M_X_double_prime.subs(s, 0)
# 分散の計算
variance = second_moment - expectation**2
# 結果を表示
display(expectation, variance)