ホーム » Qchan の投稿 (ページ 10)

作者アーカイブ: Qchan

線形変換された確率変数のモーメント母関数

投稿日:2024年4月28日　最終更新日:2024年7月28日

線形変換された確率変数のモーメント母関数を求めました。

2017 Q3(3)

投稿日:2024年4月27日　最終更新日:2024年11月4日

ポアソン分布に従う独立した2変数の和の分布を３つの方法で求めました。

コード

ポアソン分布に従う独立した2変数X1,X2と、Y=X1+X2の分布をシミュレーションしました。

# 2017 Q3(3)  2024.11.4

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
lambda_1 = 3  # X1のポアソン分布のパラメータ
lambda_2 = 4  # X2のポアソン分布のパラメータ
sample_size = 10000  # サンプルサイズ

# X1 と X2 のサンプルを生成
X1_samples = np.random.poisson(lambda_1, sample_size)
X2_samples = np.random.poisson(lambda_2, sample_size)
Y_samples = X1_samples + X2_samples  # Y = X1 + X2 のサンプル

# ヒストグラムをプロット (塗りあり)
plt.figure(figsize=(10, 4.2))  # 高さを70%に縮小
plt.hist(X1_samples, bins=range(0, 20), density=True, alpha=0.5, label=f"X1 ~ ポアソン(λ1 = {lambda_1})", color="blue")
plt.hist(X2_samples, bins=range(0, 20), density=True, alpha=0.5, label=f"X2 ~ ポアソン(λ2 = {lambda_2})", color="green")
plt.hist(Y_samples, bins=range(0, 20), density=True, alpha=0.5, label=f"Y = X1 + X2 ~ ポアソン(λ1 + λ2 = {lambda_1 + lambda_2})", color="red", histtype='stepfilled')

# ヒストグラムのカスタマイズ
plt.xlabel("値")
plt.ylabel("確率密度")
plt.title("ポアソン分布の重ね合わせ (X1, X2, および Y = X1 + X2)")
plt.legend()
plt.grid(True)
plt.show()

# CDFの計算用にPMFを求める
x_values = range(0, 20)
pmf_X1 = [np.exp(-lambda_1) * lambda_1**x / np.math.factorial(x) for x in x_values]
pmf_X2 = [np.exp(-lambda_2) * lambda_2**x / np.math.factorial(x) for x in x_values]
pmf_Y = [np.exp(-(lambda_1 + lambda_2)) * (lambda_1 + lambda_2)**x / np.math.factorial(x) for x in x_values]

# CDFを計算
X1_cdf = np.cumsum(pmf_X1)
X2_cdf = np.cumsum(pmf_X2)
Y_cdf = np.cumsum(pmf_Y)

# CDFのプロット
plt.figure(figsize=(10, 4.2))
plt.step(x_values, X1_cdf, where='mid', label=f"X1 ~ ポアソン(λ1 = {lambda_1}) CDF", color="blue")
plt.step(x_values, X2_cdf, where='mid', label=f"X2 ~ ポアソン(λ2 = {lambda_2}) CDF", color="green")
plt.step(x_values, Y_cdf, where='mid', label=f"Y = X1 + X2 ~ ポアソン(λ1 + λ2 = {lambda_1 + lambda_2}) CDF", linestyle="--", color="red")

# CDFプロットのカスタマイズ
plt.xlabel("値")
plt.ylabel("累積分布関数 (CDF)")
plt.title("ポアソン分布の累積分布関数 (X1, X2, および Y = X1 + X2)")
plt.legend()
plt.grid(True)
plt.show()

X1～Po(λ1),X1～Po(λ2)のときY=X1+X2～Po(λ1+λ2)になることをグラフの形状からも確認できました。

2017 Q3(2)

投稿日:2024年4月26日　最終更新日:2024年11月4日

ポアソン分布のモーメント母関数を求めて、それを使って期待値と分散を求めました。

コード

ポアソン分布のパラメータλを変化させシミュレーションし期待と分散が共にλになるか確認しました。

# 2017 Q3(2)  2024.11.3

import numpy as np
import matplotlib.pyplot as plt

# λの範囲を定義
lambda_values = np.arange(1, 21)
sample_size = 10000  # 各λに対するサンプル数

# サンプル平均と分散を格納するリスト
sample_means = []
sample_variances = []

# 各λに対するシミュレーションを実行
for lambda_val in lambda_values:
    # ポアソン分布のサンプルを生成
    samples = np.random.poisson(lambda_val, sample_size)
    
    # サンプルの平均と分散を計算
    sample_means.append(np.mean(samples))
    sample_variances.append(np.var(samples))

# 結果をプロット
plt.figure(figsize=(12, 6))

# サンプル平均と理論値をプロット
plt.plot(lambda_values, sample_means, label="サンプル平均 (期待値)", marker='o')
plt.plot(lambda_values, lambda_values, label="理論値 (期待値)", linestyle='--')

# サンプル分散と理論値をプロット
plt.plot(lambda_values, sample_variances, label="サンプル分散", marker='x')
plt.plot(lambda_values, lambda_values, label="理論値 (分散)", linestyle=':')

# グラフのカスタマイズ
plt.xlabel("λ")
plt.ylabel("期待値 / 分散")
plt.title("ポアソン分布の期待値と分散の可視化 (λの関数として)")
plt.legend()
plt.grid(True)
plt.show()

ポアソン分布の期待と分散が共にλになることが確認できました。

2017 Q3(1)

投稿日:2024年4月25日　最終更新日:2024年11月2日

二項分布の極限をとってポアソン分布を導出しました。

λ=5に固定しnを徐々に増加させて二項分布がポアソン分布に近づくのか確認してみます。

# 2017 Q3(1)  2024.11.2

import numpy as np
import matplotlib.pyplot as plt

# パラメータ
lambda_value = 5  # ポアソン分布の λ
n_values = [10, 50, 100, 500, 1000]  # 二項分布の異なる n 値
sample_size = 10000  # サンプル数

# サブプロットの設定
fig, axes = plt.subplots(len(n_values), 1, figsize=(8, len(n_values) * 3))
fig.suptitle("二項分布の極限がポアソン分布に収束")

# 各 n に対するシミュレーションとプロット
for i, n in enumerate(n_values):
    p = lambda_value / n  # 与えられた λ に対する p の計算
    binomial_samples = np.random.binomial(n, p, sample_size)  # 二項分布のサンプル
    poisson_samples = np.random.poisson(lambda_value, sample_size)  # ポアソン分布のサンプル（比較用）
    
    # ヒストグラムのプロット
    axes[i].hist(binomial_samples, bins=range(0, 20), alpha=0.5, label=f"二項分布(n={n}, p={p:.3f})", density=True)
    axes[i].hist(poisson_samples, bins=range(0, 20), alpha=0.5, label=f"ポアソン分布(λ={lambda_value})", density=True)
    axes[i].legend()
    axes[i].set_xlabel("x")
    axes[i].set_ylabel("確率密度")
    axes[i].set_title(f"n = {n}, p = {p:.5f}")

plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()

nが大きくなるにつれて二項分布はポアソン分布に近づくのが確認できました。

2017 Q2(4)

投稿日:2024年4月24日　最終更新日:2024年11月1日

二つの不偏推定量の分散の大きさを比較する事でどちらの推定量が望ましいかを調べました。

コード

θ=10,n=20としてシミュレーションを行い、θ’とθ’’の分布を重ねて描画してみます。

# 2017 Q2(4)  2024.11.1

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
theta_true = 10  # 真の θ の値
n = 20           # サンプルサイズ
num_trials = 1000  # シミュレーションの試行回数

# θ' = 2 * X̄ と θ'' = (n + 1) / n * X_max の推定値を記録するリスト
theta_prime_estimates = []
theta_double_prime_estimates = []

# シミュレーションを実行
for _ in range(num_trials):
    # 一様分布 U(0, theta_true) から n 個のサンプルを生成
    samples = np.random.uniform(0, theta_true, n)
    
    # θ' = 2 * X̄ を計算
    theta_prime = 2 * np.mean(samples)
    theta_prime_estimates.append(theta_prime)
    
    # θ'' = (n + 1) / n * X_max を計算
    theta_double_prime = (n + 1) / n * np.max(samples)
    theta_double_prime_estimates.append(theta_double_prime)

# ヒストグラムの表示
plt.hist(theta_prime_estimates, bins=30, edgecolor='black', density=True, alpha=0.5, label=r'$\theta\' = 2 \bar{X}$')
plt.hist(theta_double_prime_estimates, bins=30, edgecolor='black', density=True, alpha=0.5, label=r'$\theta\'\' = \frac{n + 1}{n} X_{\max}$')
plt.axvline(theta_true, color='red', linestyle='dashed', linewidth=1, label=f"真の θ = {theta_true}")
plt.xlabel('推定量')
plt.ylabel('密度')
plt.title(r'$\theta\' = 2 \bar{X}$ と $\theta\'\' = \frac{n + 1}{n} X_{\max}$ の分布')
plt.legend()
plt.show()

θ’とθ’’は同じ期待(θ)を持つものの横の広がり方が異なりθ’’の分散はθ’の分散よりも小さいことが確認できました。

次にサンプルサイズnを変化させて不偏推定量θ’とθ’’を重ねて描画してみます。

# 2017 Q2(4)  2024.11.1

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
theta_true = 10  # 真の θ の値
max_n = 100      # 最大サンプルサイズ
num_trials_per_n = 100  # 各サンプルサイズでの試行回数

# 各サンプルサイズにおける θ' と θ'' の平均を記録
theta_prime_means = []
theta_double_prime_means = []

# サンプルサイズ n を 1 から max_n まで増やしながらシミュレーション
for n in range(1, max_n + 1):
    theta_prime_estimates = []
    theta_double_prime_estimates = []
    
    for _ in range(num_trials_per_n):
        # 一様分布 U(0, theta_true) から n 個のサンプルを生成
        samples = np.random.uniform(0, theta_true, n)
        
        # θ' = 2 * X̄ を計算し、その推定値を記録
        theta_prime = 2 * np.mean(samples)
        theta_prime_estimates.append(theta_prime)
        
        # θ'' = (n + 1) / n * X_max を計算し、その推定値を記録
        theta_double_prime = (n + 1) / n * np.max(samples)
        theta_double_prime_estimates.append(theta_double_prime)
    
    # 各 n に対する θ' と θ'' の平均を保存
    theta_prime_means.append(np.mean(theta_prime_estimates))
    theta_double_prime_means.append(np.mean(theta_double_prime_estimates))

# グラフ描画
plt.plot(range(1, max_n + 1), theta_prime_means, label=r'$\mathbb{E}[\theta\']$', color='orange', alpha=0.7)
plt.plot(range(1, max_n + 1), theta_double_prime_means, label=r'$\mathbb{E}[\theta\'\']$', color='saddlebrown', alpha=0.7)
plt.axhline(theta_true, color='red', linestyle='dashed', linewidth=1, label=f"真の θ = {theta_true}")
plt.xlabel(r'サンプルサイズ $n$')
plt.ylabel('推定量の平均')
plt.title(r'サンプルサイズ $n$ と 推定量 $\theta\' = 2 \bar{X}$, $\theta\'\' = \frac{n + 1}{n} X_{\max}$ の関係')
plt.legend()
plt.show()

サンプルサイズnが増加するにつれて不偏推定量θ’とθ’’は共に真のθに近づくものの振幅が異なりθ’’の分散はθ’の分散よりも小さいことが確認できました。また収束する速度もθ’’はθ’よりも速いことが確認できました。

2017 Q2(3)

投稿日:2024年4月23日　最終更新日:2024年10月31日

サンプルデータの最大値の（n+1)/n倍が一様分布の上限値の不偏推定量であることの証明をしました。

コード

θ=10,n=20としてシミュレーションを行い、θ’’の分布を見てみます。

# 2017 Q2(3)  2024.10.31

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
theta_true = 10  # 真の θ の値
n = 20           # サンプルサイズ
num_trials = 1000  # シミュレーションの試行回数

# シミュレーションを実行
theta_double_prime_estimates = []
for _ in range(num_trials):
    # 一様分布 U(0, theta_true) から n 個のサンプルを生成
    samples = np.random.uniform(0, theta_true, n)
    # サンプルの最大値 X_max を計算し、それを用いて θ'' を計算
    theta_double_prime = (n + 1) / n * np.max(samples)
    theta_double_prime_estimates.append(theta_double_prime)

# θ'' の分布をヒストグラムで表示
plt.hist(theta_double_prime_estimates, bins=30, edgecolor='black', density=True)
plt.axvline(theta_true, color='red', linestyle='dashed', linewidth=1, label=f"真の θ = {theta_true}")
plt.xlabel(r'$\theta\'\' = \frac{n + 1}{n} X_{\max}$')
plt.ylabel('密度')
plt.title(r'$\theta\'\' = \frac{n + 1}{n} X_{\max}$ の分布')
plt.legend()
plt.show()

θ’’は、元となるXmaxの分布を右にずらした形状を取っています。期待値は真のθに一致しているようです。

次にサンプルサイズnを変化させて不偏推定量θ’’がどうなるのか確認をします。

# 2017 Q2(3)  2024.10.31

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
theta_true = 10  # 真の θ の値
max_n = 100      # 最大サンプルサイズ
num_trials_per_n = 100  # 各サンプルサイズでの試行回数

# 各サンプルサイズにおける θ'' の平均を記録
theta_double_prime_means = []

# サンプルサイズ n を 1 から max_n まで増やしながらシミュレーション
for n in range(1, max_n + 1):
    theta_double_prime_estimates = []
    for _ in range(num_trials_per_n):
        # 一様分布 U(0, theta_true) から n 個のサンプルを生成
        samples = np.random.uniform(0, theta_true, n)
        # サンプルの最大値 X_max を計算し、それを用いて θ'' を計算
        theta_double_prime = (n + 1) / n * np.max(samples)
        theta_double_prime_estimates.append(theta_double_prime)
    # 各 n に対する θ'' の平均を保存
    theta_double_prime_means.append(np.mean(theta_double_prime_estimates))

# グラフ描画
plt.plot(range(1, max_n + 1), theta_double_prime_means, label=r'$\mathbb{E}[\theta\'\']$')
plt.axhline(theta_true, color='red', linestyle='dashed', linewidth=1, label=f"真の θ = {theta_true}")
plt.xlabel(r'サンプルサイズ $n$')
plt.ylabel(r'推定量 $\theta\'\' = \frac{n + 1}{n} X_{\max}$ の平均')
plt.title(r'サンプルサイズ $n$ と $\theta\'\' = \frac{n + 1}{n} X_{\max}$ の関係')
plt.legend()
plt.show()

サンプルサイズnが増加するにつれて不偏推定量θ’’は真のθに近づくことが確認できました。

2017 Q2(2)

投稿日:2024年4月22日　最終更新日:2024年10月30日

標本平均の２倍が一様分布の上限値の不偏推定量であることの証明をしました。

コード

θ=10,n=20としてシミュレーションを行い、θ’の分布を見てみます。

# 2017 Q2(2)  2024.10.30

import numpy as np
import matplotlib.pyplot as plt

# シミュレーションのパラメータ
theta_true = 10  # 真の θ の値
n = 20           # 1試行あたりのサンプル数
num_trials = 1000  # シミュレーションの試行回数

# シミュレーションを実行
theta_prime_estimates = []
for _ in range(num_trials):
    # 一様分布 U(0, theta_true) から n 個のサンプルを生成
    samples = np.random.uniform(0, theta_true, n)
    # サンプル平均 X̄ を計算し、それを用いて θ' を計算
    theta_prime = 2 * np.mean(samples)
    theta_prime_estimates.append(theta_prime)

# θ' の分布をヒストグラムで表示
plt.hist(theta_prime_estimates, bins=30, edgecolor='black', density=True)
plt.axvline(theta_true, color='red', linestyle='dashed', linewidth=1, label=f"真の θ = {theta_true}")
plt.xlabel(r'$\theta\' = 2\bar{X}$')
plt.ylabel('密度')
plt.title(r'$\theta\' = 2\bar{X}$ の分布')
plt.legend()
plt.show()

θ’は真のθを中心に左右対称にバラついています。θ’は不偏であるように見えます。

次にサンプルサイズnを変化させて不偏推定量θ’がどうなるのか確認をします。

# 2017 Q2(2)  2024.10.30

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
theta_true = 10  # 真の θ の値
max_n = 100      # 最大サンプルサイズ
num_trials_per_n = 100  # 各サンプルサイズでの試行回数

# 各サンプルサイズにおける θ' の平均を記録
theta_prime_means = []

# サンプルサイズ n を 1 から max_n まで増やしながらシミュレーション
for n in range(1, max_n + 1):
    theta_prime_estimates = []
    for _ in range(num_trials_per_n):
        # 一様分布 U(0, theta_true) から n 個のサンプルを生成
        samples = np.random.uniform(0, theta_true, n)
        # サンプル平均 X̄ を計算し、それを用いて θ' を計算
        theta_prime = 2 * np.mean(samples)
        theta_prime_estimates.append(theta_prime)
    # 各 n に対する θ' の平均を保存
    theta_prime_means.append(np.mean(theta_prime_estimates))

# グラフ描画
plt.plot(range(1, max_n + 1), theta_prime_means, label=r'$\mathbb{E}[\theta\']$')
plt.axhline(theta_true, color='red', linestyle='dashed', linewidth=1, label=f"真の θ = {theta_true}")
plt.xlabel(r'サンプルサイズ $n$')
plt.ylabel(r'推定量 $\theta\' = 2 \bar{X}$ の平均')
plt.title(r'サンプルサイズ $n$ と $\theta\' = 2 \bar{X}$ の関係')
plt.legend()
plt.show()

サンプルサイズnが増加するにつれて不偏推定量θ’は真のθに近づくことが確認できました。

2017 Q2(1)

投稿日:2024年4月21日　最終更新日:2024年10月30日

一様分布の上限の最尤推定量を求めました。

コード

θ=10,n=20としてシミュレーションを行い、 $\hat{\theta}$ の分布を見てみます。

# 2017 Q2(1)  2024.10.29

import numpy as np
import matplotlib.pyplot as plt

# シミュレーションのパラメータ
theta_true = 10  # 真の θ の値
n = 20           # 1試行あたりのサンプル数
num_trials = 1000  # シミュレーションの試行回数

# シミュレーションを実行
theta_estimates = []
for _ in range(num_trials):
    # 一様分布 U(0, theta_true) から n 個のサンプルを生成
    samples = np.random.uniform(0, theta_true, n)
    # サンプルの最大値を推定値 θ_hat として記録
    theta_hat = np.max(samples)
    theta_estimates.append(theta_hat)

# 推定値の分布をヒストグラムで表示
plt.hist(theta_estimates, bins=30, edgecolor='black', density=True)
plt.axvline(theta_true, color='red', linestyle='dashed', linewidth=1, label=f"真の θ = {theta_true}")
plt.xlabel(r'$\hat{\theta}$')
plt.ylabel('密度')
plt.title(r'$\hat{\theta} = X_{\max}$ の分布')
plt.legend()
plt.show()

$\hat{\theta}$ は真のθより大きくならないため、結果としてθよりやや小さな値を取っているように見えます。

次にサンプルサイズnを変化させて最尤推定量 $\hat{\theta}$ がどうなるのか確認をします。

# 2017 Q2(1)  2024.10.29

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
theta_true = 10  # 真の θ の値
max_n = 100  # 最大サンプルサイズ
num_trials_per_n = 100  # 各サンプルサイズでの試行回数

# 各サンプルサイズにおける θ^ の平均を記録
theta_hat_means = []

# サンプルサイズ n を 1 から max_n まで増やしながらシミュレーション
for n in range(1, max_n + 1):
    theta_estimates = []
    for _ in range(num_trials_per_n):
        # 一様分布 U(0, theta_true) から n 個のサンプルを生成
        samples = np.random.uniform(0, theta_true, n)
        # サンプルの最大値を推定値 θ_hat として記録
        theta_hat = np.max(samples)
        theta_estimates.append(theta_hat)
    # 各 n に対する θ^ の平均を保存
    theta_hat_means.append(np.mean(theta_estimates))

# グラフ描画
plt.plot(range(1, max_n + 1), theta_hat_means, label=r'$\mathbb{E}[\hat{\theta}]$')
plt.axhline(theta_true, color='red', linestyle='dashed', linewidth=1, label=f"真の θ = {theta_true}")
plt.xlabel(r'サンプルサイズ $n$')
plt.ylabel(r'推定量 $\hat{\theta}$ の平均')
plt.title(r'サンプルサイズ $n$ と $\hat{\theta}$ の関係')
plt.legend()
plt.show()

サンプルサイズnが増加するにつれて最尤推定量 $\hat{\theta}$ は真のθに近づくことが確認できました。

2017 Q1(5)

投稿日:2024年4月20日　最終更新日:2024年10月28日

正規分布の分散の最尤推定量を母平均が既知の場合と未知の場合でそれぞれ求めましまた。

コード

サンプルサイズnを50に固定し尤度関数を可視化します。また最尤推定量 $T^2 = \frac{1}{n} \sum_{i=1}^n (X_i - \mu)^2$ と $S^2 = \frac{1}{n-1} \sum_{i=1}^n (X_i - \bar{X})^2$ をプロットしてみます。

# 2017 Q1(5)  2024.10.28

# 最尤推定量 (T^2, S^2) の計算と尤度関数の可視化
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# シミュレーションパラメータの設定
mu = 5.0             # 母平均（既知）
sigma = 2.0          # 母標準偏差（√母分散）
selected_sample_size = 50  # 固定サンプルサイズ
true_variance = sigma**2    # 母分散

# 正規分布からのサンプル生成（サンプルサイズは selected_sample_size で固定）
sample = np.random.normal(mu, sigma, selected_sample_size)  # 固定サンプル生成

# 最尤推定量の計算
T2 = np.sum((sample - mu) ** 2) / selected_sample_size  # 最尤推定量 T^2 の計算 (μが既知)
sample_mean = np.mean(sample)
S2 = np.sum((sample - sample_mean) ** 2) / (selected_sample_size - 1)  # 最尤推定量 S^2 の計算 (μが未知)

# 尤度関数のプロットのためのパラメータ設定
sigma_squared_trial_values = np.linspace(0.1, 5.0, 100)  # σ^2 の値を試行的に設定
likelihoods = []  # 尤度関数の値を保存するリスト

# 尤度関数の計算（固定サンプルに基づく）
for sigma_squared_trial in sigma_squared_trial_values:
    likelihood = np.prod(norm.pdf(sample, loc=mu, scale=np.sqrt(sigma_squared_trial)))  # 尤度関数の値を計算
    likelihoods.append(likelihood)

# グラフの作成
plt.figure(figsize=(10, 6))
plt.plot(sigma_squared_trial_values, likelihoods, label='尤度関数', color='blue')
plt.axvline(x=T2, color='green', linestyle='--', label='$T^2$ (最尤推定量, μが既知)', linewidth=2)
plt.axvline(x=S2, color='orange', linestyle='--', label='$S^2$ (最尤推定量, μが未知)', linewidth=2)
plt.axvline(x=true_variance, color='red', linestyle='-.', label='母分散 $\\sigma^2$', linewidth=2)
plt.title(f'サンプルサイズ {selected_sample_size} における尤度関数と推定量の比較')
plt.xlabel('$\sigma^2$ の試行値')
plt.ylabel('尤度')
plt.legend()
plt.grid(True)

# グラフの表示
plt.tight_layout()
plt.show()

最尤推定量 $T^2$ と $S^2$ は尤度が最大となる $\sigma^2$ の値を取っています。

もう一度実行してみます。

サンプルによっては最尤推定量 $T^2$ と $S^2$ は $\sigma^2$ から乖離する場合もあります。

次に、サンプルサイズを大きくすることで母分散に近づくことを確認します。

# 2017 Q1(5)  2024.10.28

from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt

# シミュレーションパラメータの設定
mu = 5.0             # 母平均（既知）
sigma = 2.0          # 母標準偏差（√母分散）
true_variance = sigma**2    # 母分散
sample_sizes = range(50, 10050, 50)  # サンプルサイズを設定

# シミュレーション結果を保存するリスト
T2_estimates = []  # μが既知の場合の推定量 T^2
S2_estimates = []  # μが未知の場合の推定量 S^2

# 各サンプルサイズでのシミュレーションの実行
for n in sample_sizes:
    # 正規分布からサンプルを抽出
    sample = np.random.normal(mu, sigma, n)
    
    # μが既知の場合の最尤推定量 T^2
    T2 = np.sum((sample - mu) ** 2) / n
    
    # μが未知の場合の最尤推定量 S^2
    sample_mean = np.mean(sample)
    S2 = np.sum((sample - sample_mean) ** 2) / (n - 1)
    
    # 推定量をリストに保存
    T2_estimates.append(T2)
    S2_estimates.append(S2)

# グラフの作成（理論値を追加）
plt.figure(figsize=(12, 8))

# μが既知の場合のプロット (緑)
plt.plot(sample_sizes, T2_estimates, marker='o', color='green', label='$T^2$ (μが既知の場合)', linestyle='--')

# μが未知の場合のプロット (オレンジ)
plt.plot(sample_sizes, S2_estimates, marker='x', color='orange', label='$S^2$ (μが未知の場合)', linestyle='--')

# 理論値の線を追加 (赤)
plt.axhline(y=true_variance, color='red', linestyle='-.', label='理論値: $\sigma^2$', linewidth=2)

# グラフの設定
plt.title('最尤推定量 $T^2$ と $S^2$ のサンプルサイズごとの結果')
plt.xlabel('サンプルサイズ $n$')
plt.ylabel('推定量')
plt.legend()
plt.grid(True)

# グラフの表示
plt.tight_layout()
plt.show()

サンプルサイズを大きくすることで最尤推定量 $T^2$ と $S^2$ は母分散 $\sigma^2$ に近づきました。

2017 Q1(4)

投稿日:2024年4月19日　最終更新日:2024年10月27日

サンプル数が十分に大きい時のサンプル平均の歪度と尖度がどのようになるか考察した。

コード

Xの標本平均の歪度と尖度（過剰尖度）はサンプルサイズnが十分大きいときにどう変化するのか確認します。

# 2017 Q1(4)  2024.10.27

from scipy.stats import skew, kurtosis, gamma
import numpy as np
import matplotlib.pyplot as plt

# シミュレーションパラメータの設定
#np.random.seed(42)  # 再現性のためのシード
num_simulations = 300000  # シミュレーションの繰り返し回数
sample_sizes = range(10, 1100, 50)  # サンプルサイズを広げて設定

# 歪度と尖度を持つ母集団分布としてガンマ分布を使用（shape=2で正の歪度と尖度を持つ）
shape_param = 2.0  # ガンマ分布のshapeパラメータ
scale_param = 2.0  # ガンマ分布のscaleパラメータ
population_dist = gamma(a=shape_param, scale=scale_param)

# シミュレーション結果を保存するリスト
sample_mean_skewness = []
sample_mean_kurtosis = []

# 各サンプルサイズでのシミュレーションの実行
for n in sample_sizes:
    skew_values = []
    kurt_values = []
    
    # シミュレーションを繰り返し
    for _ in range(num_simulations):
        sample = population_dist.rvs(size=n)  # サンプルを抽出
        sample_mean = np.mean(sample)  # 標本平均
        
        # 標本平均をリストに保存
        skew_values.append(sample_mean)
        kurt_values.append(sample_mean)
    
    # 標本平均の歪度と尖度を計算
    skewness_of_sample_mean = skew(skew_values)
    kurtosis_of_sample_mean = kurtosis(kurt_values, fisher=False) - 3  # 過剰尖度に変換
    
    # 結果を保存
    sample_mean_skewness.append(skewness_of_sample_mean)
    sample_mean_kurtosis.append(kurtosis_of_sample_mean)

# グラフの作成
plt.figure(figsize=(12, 6))

# 歪度のプロット
plt.subplot(1, 2, 1)
plt.plot(sample_sizes, sample_mean_skewness, marker='o', label='歪度')
plt.axhline(y=0, color='r', linestyle='--', label='理論値: 0')
plt.title('サンプルサイズと歪度の関係')
plt.xlabel('サンプルサイズ $n$')
plt.ylabel('歪度 ($\\bar{X}$)')
plt.legend()
plt.grid(True)

# 尖度のプロット
plt.subplot(1, 2, 2)
plt.plot(sample_sizes, sample_mean_kurtosis, marker='o', label='尖度（過剰尖度）')
plt.axhline(y=0, color='r', linestyle='--', label='理論値: 0')
plt.title('サンプルサイズと尖度の関係')
plt.xlabel('サンプルサイズ $n$')
plt.ylabel('尖度 ($\\bar{X}$)')
plt.legend()
plt.grid(True)

# グラフの表示
plt.tight_layout()
plt.show()

Xの標本平均の尖度（過剰尖度）はサンプルサイズnが大きくなると0に近づくことが確認できました。しかし歪度は0に近づくスピードは遅いです。Xが従うガンマ分布は元々歪度を持つ非対称な分布なので、そのようになります。サンプルサイズnを10,000まで大きくしてみます。

nが十分に大きくなったので歪度も0に近づきました。

次に、サンプルサイズnに伴う標本平均の分布が正規分布に近づく様子を確認します。

# 2017 Q1(4)  2024.10.27

import seaborn as sns
from scipy.stats import norm

# サンプルサイズの設定
selected_sample_sizes = [5, 20, 100, 200]  # 表示するサンプルサイズの選択
num_simulations = 10000  # シミュレーションの繰り返し回数

# グラフの作成
plt.figure(figsize=(14, 10))

# 各サンプルサイズでのシミュレーションとヒストグラムの描画
for idx, n in enumerate(selected_sample_sizes):
    sample_means = []
    
    # シミュレーションの実行
    for _ in range(num_simulations):
        sample = population_dist.rvs(size=n)  # サンプルを抽出
        sample_mean = np.mean(sample)  # 標本平均を計算
        sample_means.append(sample_mean)
    
    # ヒストグラムのプロット
    plt.subplot(2, 2, idx + 1)
    sns.histplot(sample_means, bins=30, kde=True, stat="density", color='blue', label='標本平均の分布')
    
    # 正規分布のプロット
    mean_of_means = np.mean(sample_means)
    std_of_means = np.std(sample_means)
    x = np.linspace(min(sample_means), max(sample_means), 100)
    plt.plot(x, norm.pdf(x, mean_of_means, std_of_means), 'r--', label='正規分布')
    
    # グラフの設定
    plt.title(f'サンプルサイズ {n} の標本平均の分布')
    plt.xlabel('標本平均')
    plt.ylabel('密度')
    plt.legend()
    plt.grid(True)

# グラフの表示
plt.tight_layout()
plt.show()

サンプルサイズnが大きくなるほど標本平均の分布は正規分布に近づくことが確認できました。

作者アーカイブ: Qchan

カテゴリー

投稿一覧