ホーム » 不偏推定量
「不偏推定量」カテゴリーアーカイブ
2017 Q2(4)
二つの不偏推定量の分散の大きさを比較する事でどちらの推定量が望ましいかを調べました。
コード
θ=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)
サンプルデータの最大値の（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)
標本平均の２倍が一様分布の上限値の不偏推定量であることの証明をしました。
コード
θ=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が増加するにつれて不偏推定量θ’は真のθに近づくことが確認できました。
2018 Q1(1)
不偏分散が母分散の不偏推定量であることを示しました。
コード
nを2～100に変化させて、不偏分散と標本分散を比較してみます。
import numpy as np
import matplotlib.pyplot as plt
# パラメータ設定
mu = 0 # 母集団平均
sigma = 2 # 母集団標準偏差
sigma_squared = sigma ** 2 # 真の母分散
n_values = range(2, 101, 2) # サンプルサイズ n を2から100まで2ステップで変化させる
num_trials = 3000 # 各 n に対して100回の試行を行う
# 不偏分散 (1/(n-1)) と 標本分散 (1/n) を計算するためのリスト
unbiased_variances = []
biased_variances = []
# 各サンプルサイズ n で分散を計算
for n in n_values:
unbiased_variance_sum = 0
biased_variance_sum = 0
# 各サンプルサイズ n に対して複数回試行して平均を計算
for _ in range(num_trials):
# 正規分布に従うサンプルを生成
sample = np.random.normal(mu, sigma, n)
# 標本平均
sample_mean = np.mean(sample)
# 不偏分散 (1/(n-1))
unbiased_variance = np.sum((sample - sample_mean) ** 2) / (n - 1)
unbiased_variance_sum += unbiased_variance
# 標本分散 (1/n)
biased_variance = np.sum((sample - sample_mean) ** 2) / n
biased_variance_sum += biased_variance
# 各 n に対する平均分散をリストに追加
unbiased_variances.append(unbiased_variance_sum / num_trials)
biased_variances.append(biased_variance_sum / num_trials)
# グラフを描画
plt.plot(n_values, unbiased_variances, label="不偏分散 (1/(n-1))", color='blue', marker='o')
plt.plot(n_values, biased_variances, label="標本分散 (1/n)", color='red', linestyle='--', marker='x')
# 真の分散を水平線で描画
plt.axhline(y=sigma_squared, color='green', linestyle='-', label=f'真の分散 = {sigma_squared}')
# グラフの設定
plt.title('サンプルサイズに対する標本分散と不偏分散の比較')
plt.xlabel('サンプルサイズ n')
plt.ylabel('分散')
plt.legend()
plt.grid(True)
plt.show()
標本分散には不偏性はなく、不偏分散には不偏性があることが確認できました。
2019 Q3(4)
一様分布に従う複数の確率変数がとる最大値の期待値と、その最大値から上限θの不偏推定量を求めました。
コード
数式を使って期待値E(Y)とθの不偏推定量を求めます。
# 2019 Q3(4) 2024.9.20
import sympy as sp
# シンボリック変数の定義
y, theta, n = sp.symbols('y theta n', positive=True, real=True)
# 最大値Yの確率密度関数 f_Y(y)
f_Y = (n / theta**n) * y**(n-1)
# 期待値 E[Y] の定義
E_Y = sp.integrate(y * f_Y, (y, 0, theta))
# 期待値 E[Y] の表示
print(f"E[Y] の導出結果:")
display(sp.simplify(E_Y))
# 改行を追加
print()
# 新しいシンボリック変数Yの定義
Y = sp.symbols('Y', positive=True, real=True)
# 理論的な E[Y] の式（E[Y] = n/(n+1) * theta）
theoretical_E_Y = (n / (n + 1)) * theta
# Y を基にして theta を求める方程式を解く
theta_hat = sp.solve(Y - theoretical_E_Y, theta)[0]
# 不偏推定量の結果表示
print(f"不偏推定量 θ̂ の導出結果:")
display(theta_hat)
式の形は少し違いますが正しいですね。
次に、数値シミュレーションにより期待値E(Y)とθの不偏推定量求め理論値と比較します。
# 2019 Q3(4) 2024.9.20
import numpy as np
import matplotlib.pyplot as plt
# パラメータ設定
theta_true = 10 # 真のθ
n = 15 # サンプルサイズ
num_simulations = 10000 # シミュレーション回数
# シミュレーションでの最大値Yを格納するリスト
max_values = []
# シミュレーション開始
for _ in range(num_simulations):
# U(0, theta_true) から n 個の乱数を生成
samples = np.random.uniform(0, theta_true, n)
# その中での最大値Yを記録
max_values.append(np.max(samples))
# シミュレーションによる期待値E[Y]を計算
simulated_E_Y = np.mean(max_values)
# 理論値E[Y]を計算
theoretical_E_Y = (n / (n + 1)) * theta_true
# 不偏推定量のシミュレーション
theta_hat_simulated = (n + 1) / n * np.mean(max_values)
# 結果の表示
print(f"シミュレーションによる E[Y]: {simulated_E_Y}")
print(f"理論値 E[Y]: {theoretical_E_Y}")
print(f"シミュレーションによる不偏推定量 θ̂: {theta_hat_simulated}")
print(f"真の θ: {theta_true}")
# ヒストグラムを描画して結果を視覚化
plt.hist(max_values, bins=30, density=True, alpha=0.7, color='blue', edgecolor='black', label='シミュレーション結果')
plt.axvline(simulated_E_Y, color='red', linestyle='--', label=f'シミュレーション E[Y] = {simulated_E_Y:.2f}')
plt.axvline(theoretical_E_Y, color='green', linestyle='--', label=f'理論値 E[Y] = {theoretical_E_Y:.2f}')
# グラフ設定
plt.title('最大値Yの分布とE[Y]')
plt.xlabel('Y の値')
plt.ylabel('密度')
plt.legend()
plt.grid(True)
plt.show()シミュレーションによる E[Y]: 9.384118720046189
理論値 E[Y]: 9.375
シミュレーションによる不偏推定量 θ̂: 10.009726634715935
真の θ: 10
期待値E(Y)とθの不偏推定量は理論値に近い値を取りました。