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

作者アーカイブ: Qchan

投稿一覧

2019 Q5(1)-1

ラプラス分布の期待値を求めました。

 

コード

ξ=0,λ=1のラプラス分布 g(\mu) = \frac{\lambda}{2} \exp\left(-\lambda |\mu - \xi|\right) の期待値E[μ]と分散V[μ]を求めます。

# 2019 Q5(1)  2024.9.27

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
xi = 0  # 理論的な期待値
lambda_ = 1  # スケールパラメータ
n_samples = 10000  # サンプル数

# ラプラス分布からランダムサンプルを生成
samples = np.random.laplace(loc=xi, scale=1/lambda_, size=n_samples)

# シミュレーションの期待値と分散を計算
sample_mean = np.mean(samples)
sample_variance = np.var(samples)

# 理論値
theoretical_mean = xi
theoretical_variance = 2 / lambda_**2

# 理論値とシミュレーション値を出力
print(f"シミュレーションによる期待値 E[μ]: {sample_mean}")
print(f"理論的な期待値 E[μ]: {theoretical_mean}")
print(f"シミュレーションによる分散 V[μ]: {sample_variance}")
print(f"理論的な分散 V[μ]: {theoretical_variance}")

# 理論的なラプラス分布のPDFを描画
mu_values = np.linspace(-10, 10, 1000)
pdf_values = (lambda_ / 2) * np.exp(-lambda_ * np.abs(mu_values - xi))

# ヒストグラムの描画
plt.hist(samples, bins=50, density=True, alpha=0.6, color='g', label='シミュレーションデータ')

# 理論的なPDFを線で描画
plt.plot(mu_values, pdf_values, 'r-', lw=2, label='理論的なPDF')

# グラフの設定
plt.title('ラプラス分布:シミュレーション vs 理論的なPDF')
plt.xlabel('値')
plt.ylabel('密度')
plt.legend()
plt.show()
シミュレーションによる期待値 E[μ]: -8.093411765355611e-05
理論的な期待値 E[μ]: 0
シミュレーションによる分散 V[μ]: 2.043434597873531
理論的な分散 V[μ]: 2.0

シミュレーションによる期待値E[μ]と分散V[μ]は理論値に近い値をとりました。

2019 Q4(4)

コーシー分布の検定で、与えられた棄却域が最強力検定になることを示しました。

 

コード

棄却域をR={x:1<x<3}にする検定が最強力検定になるのか確認するために、尤度比λ(x)のグラフを描画します。

# 2019 Q4(4)  2024.9.26

import numpy as np
import matplotlib.pyplot as plt

# 尤度比 λ(x) の計算
def likelihood_ratio(x):
    return (1 + x**2) / (1 + (x - 1)**2)

# x の範囲を設定
x_values = np.linspace(-5, 5, 500)

# λ(x) の値を計算
lambda_values = likelihood_ratio(x_values)

# グラフを描画
plt.figure(figsize=(10, 6))
plt.plot(x_values, lambda_values, label=r'尤度比 $\lambda(x) = \frac{1 + x^2}{1 + (x - 1)^2}$', color='b')

# 棄却域 1 < x < 3 を塗りつぶす(透明度を追加)
plt.axvspan(1, 3, color='yellow', alpha=0.3, label=r'棄却域 $1 < x < 3$')

# λ(x) >= 2 の部分を赤で塗る
plt.fill_between(x_values, lambda_values, 2, where=(lambda_values >= 2), color='red', alpha=0.3, label=r'$\lambda(x) \geq 2$')

# グラフの装飾
plt.axhline(2, color='green', linestyle='--', label=r'閾値 $\lambda(x) = 2$')
plt.axhline(0, color='black',linewidth=0.5)
plt.axvline(0, color='black',linewidth=0.5)
plt.title(r'尤度比 $\lambda(x)$ と棄却域 $1 < x < 3$, $\lambda(x) \geq 2$ の範囲')
plt.xlabel('x')
plt.ylabel(r'$\lambda(x)$')
plt.grid(True)
plt.legend()

# グラフを表示
plt.show()

棄却域は尤度比λ(x)>c (ここではc=2)で定義できるので、これは最強力検定と言えます。

ところでx->-∞で尤度比λ(x)はどこに収束するのでしょうか。

import numpy as np
import matplotlib.pyplot as plt

# 尤度比 λ(x) の計算
def likelihood_ratio(x):
    return (1 + x**2) / (1 + (x - 1)**2)

# x の範囲を設定(左側を広げる)
x_values = np.linspace(-50, 5, 500)

# λ(x) の値を計算
lambda_values = likelihood_ratio(x_values)

# グラフを描画
plt.figure(figsize=(10, 6))
plt.plot(x_values, lambda_values, label=r'尤度比 $\lambda(x) = \frac{1 + x^2}{1 + (x - 1)^2}$', color='b')

# 閾値を視覚化
plt.axhline(1, color='purple', linestyle='--', label=r'極限 $\lambda(x) \to 1$ as $x \to -\infty$')
plt.axhline(0, color='black',linewidth=0.5)
plt.axvline(0, color='black',linewidth=0.5)
plt.title(r'尤度比 $\lambda(x)$ の概形 ($x \to -\infty$ の範囲)')
plt.xlabel('x')
plt.ylabel(r'$\lambda(x)$')
plt.grid(True)
plt.legend()

# グラフを表示
plt.show()

1.0に収束するようです。

よって、cは1以上で、且つ尤度比λ(x)が最大となるx = 1.61803398874989(黄金比)が棄却域に含まれていれば最強力検定になります。

2019 Q4(3)

コーシー分布のパラメータ違いの尤度比の概形を描きました。

 

コード

尤度比λ(x)のグラフを描画します。

# 2019 Q4(3)  2024.9.25

import numpy as np
import matplotlib.pyplot as plt

# 尤度比 λ(x) の計算
def likelihood_ratio(x):
    return (1 + x**2) / (1 + (x - 1)**2)

# x の範囲を設定
x_values = np.linspace(-5, 5, 500)

# λ(x) の値を計算
lambda_values = likelihood_ratio(x_values)

# グラフを描画
plt.figure(figsize=(10, 6))
plt.plot(x_values, lambda_values, label=r'尤度比 $\lambda(x) = \frac{1 + x^2}{1 + (x - 1)^2}$', color='b')
plt.scatter([1, 3], [2, 2], color='red', zorder=5)  # x = 1 および x = 3 の点をプロット

# グラフの装飾
plt.axhline(0, color='black',linewidth=0.5)
plt.axvline(0, color='black',linewidth=0.5)
plt.title(r'尤度比 $\lambda(x)$ の概形')
plt.xlabel('x')
plt.ylabel(r'$\lambda(x)$')
plt.grid(True)
plt.legend()

# グラフを表示
plt.show()

谷と山を持つ形状をしています。

次に、尤度関数λ(x)の1次導関数と2次導関数から、極値と形状(山か谷か)を導きます。

# 2019 Q4(3)  2024.9.25

import sympy as sp
from IPython.display import display

# 尤度比 λ(x) の式を定義
x = sp.symbols('x')
lambda_x = (1 + x**2) / (1 + (x - 1)**2)

# 1次導関数を計算
lambda_prime = sp.diff(lambda_x, x)
print("1次導関数 λ'(x):")
display(lambda_prime)
print()

# 1次導関数が0となる点(極値候補)を解く
critical_points = sp.solve(lambda_prime, x)
print("極値候補:")
display(critical_points)
print()

# 2次導関数を計算して極値の性質を確認
lambda_double_prime = sp.diff(lambda_prime, x)
print("2次導関数 λ''(x):")
display(lambda_double_prime)
print()

# 各極値候補の 2次導関数の値を調べる
for point in critical_points:
    second_derivative_value = lambda_double_prime.subs(x, point)
    
    # 極値の場所を数値評価する
    point_value = point.evalf()
    
    # 2次導関数の値も数値評価する
    second_derivative_value_numeric = second_derivative_value.evalf()
    
    if second_derivative_value > 0:
        print(f"x = {point} で谷(極小値) (数値: {point_value})")
        print(f"2次導関数の値: {second_derivative_value_numeric}")
    elif second_derivative_value < 0:
        print(f"x = {point} で山(極大値) (数値: {point_value})")
        print(f"2次導関数の値: {second_derivative_value_numeric}")
    else:
        print(f"x = {point} で特異点か平坦な極値 (数値: {point_value})")
        print(f"2次導関数の値: {second_derivative_value_numeric}")

以上のように求まりました。

2019 Q4(2)

コーシー分布の検出力を計算しました。

 

コード

数式を使って検出力1-βを求めます。

# 2019 Q4(2)  2024.9.24

import numpy as np
from scipy.integrate import quad

# コーシー分布の確率密度関数 (theta = 1)
def cauchy_pdf_theta_1(x):
    return 1 / (np.pi * (1 + (x - 1)**2))

# 積分範囲(棄却域R = (1, 3))
lower_bound = 1
upper_bound = 3

# 積分を実行
power, error = quad(cauchy_pdf_theta_1, lower_bound, upper_bound)

# 結果を表示(小数第3位まで表示)
print(f"検出力 (1 - β): {power:.3f}")
検出力 (1 - β): 0.352

手計算と一致しました。

次に、数値シミュレーションで計算をしてみます。なお、コーシー分布は裾が重いため-8 ~8の範囲になるようにフィルターを掛けることにします。

# 2019 Q4(2)  2024.9.24

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

# シミュレーションのパラメータ
np.random.seed(43)
num_trials = 100000

# 棄却域 R = (1, 3)
lower_bound = 1
upper_bound = 3

# コーシー分布 (θ = 1) からのサンプルを生成
samples = cauchy.rvs(loc=1, scale=1, size=num_trials)

# -8 から 8 の範囲にサンプルを制限して外れ値を除外
filtered_samples = samples[(samples > -8) & (samples < 8)]

# 棄却域に入っているかどうかを判定
reject = (samples > lower_bound) & (samples < upper_bound)

# 棄却域に入った割合が検出力 (1 - β)
power_simulated = np.mean(reject)

# 結果を表示
print(f"シミュレーションによる検出力 (1 - β): {power_simulated:.3f}")

# ヒストグラムのプロット(フィルタリングしたサンプルを使う)
plt.figure(figsize=(10, 6))
plt.hist(filtered_samples, bins=100, density=True, alpha=0.5, color='blue', label='シミュレーション結果(フィルタリング済み)')

# 理論的なコーシー分布の確率密度関数 (theta = 1) をプロット
x = np.linspace(-8, 8, 1000)
pdf_theta_1 = cauchy.pdf(x, loc=1)  # 対立仮説の理論曲線 (θ = 1)
plt.plot(x, pdf_theta_1, 'r-', lw=2, label='理論的コーシー分布 (θ = 1)')

# 理論的なコーシー分布の確率密度関数 (theta = 0) をプロット
pdf_theta_0 = cauchy.pdf(x, loc=0)  # 帰無仮説の理論曲線 (θ = 0)
plt.plot(x, pdf_theta_0, 'g--', lw=2, label='理論的コーシー分布 (θ = 0)')

# 棄却域を塗りつぶす
plt.axvspan(lower_bound, upper_bound, color='red', alpha=0.3, label='棄却域 (1 < x < 3)')

# グラフの設定
plt.xlim(-8, 8)
plt.xlabel('x')
plt.ylabel('密度')
plt.title(f'コーシー分布のシミュレーションと棄却域\nシミュレーションによる検出力 (1 - β) = {power_simulated:.3f}')
plt.legend()
plt.grid(True)

# グラフを表示
plt.show()
シミュレーションによる検出力 (1 - β): 0.355

シミュレーション結果も手計算とほぼ一致しました。

2019 Q4(1)

コーシー分布の検定での第一種過誤確率を求めました。

 

コード

数式を使って第一種の過誤確率αを求めます。

# 2019 Q4(1)  2024.9.23

import numpy as np
from scipy.integrate import quad

# コーシー分布の確率密度関数 (theta = 0)
def cauchy_pdf(x):
    return 1 / (np.pi * (1 + x**2))

# 積分範囲(棄却域R = (1, 3))
lower_bound = 1
upper_bound = 3

# 積分を実行
alpha, error = quad(cauchy_pdf, lower_bound, upper_bound)

# 結果を表示(小数第3位まで表示)
print(f"第一種の過誤確率 α: {alpha:.3f}")
第一種の過誤確率 α: 0.148

手計算と一致しました。

次に、数値シミュレーションで計算をしてみます。なお、コーシー分布は裾が重いため-8 ~8の範囲になるようにフィルターを掛けることにします。

# 2019 Q4(1)  2024.9.23

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

# シミュレーションのパラメータ
np.random.seed(43)
num_trials = 100000

# 棄却域 R = (1, 3)
lower_bound = 1
upper_bound = 3

# コーシー分布 (θ = 0) からのサンプルを生成
samples = np.random.standard_cauchy(size=num_trials)

# -8 から 8 の範囲にサンプルを制限して外れ値を除外
filtered_samples = samples[(samples > -8) & (samples < 8)]

# 棄却域に入っているかどうかを判定
reject = (samples > lower_bound) & (samples < upper_bound)

# 棄却域に入った割合が第一種の過誤確率 α
alpha_simulated = np.mean(reject)

# 結果を表示
print(f"シミュレーションによる第一種の過誤確率 α: {alpha_simulated:.3f}")

# ヒストグラムのプロット(フィルタリングしたサンプルを使う)
plt.figure(figsize=(10, 6))
plt.hist(filtered_samples, bins=100, density=True, alpha=0.5, color='blue', label='シミュレーション結果(フィルタリング済み)')

# 理論的なコーシー分布の確率密度関数 (theta = 0) をプロット
x = np.linspace(-8, 8, 1000)  # 横のレンジを -8 から 8 に設定
pdf = cauchy.pdf(x)
plt.plot(x, pdf, 'r-', lw=2, label='理論的コーシー分布の確率密度関数')

# 棄却域を塗りつぶす
plt.axvspan(lower_bound, upper_bound, color='red', alpha=0.3, label='棄却域 (1 < x < 3)')

# グラフの設定
plt.xlim(-8, 8)  # 横のレンジを -8 から 8 に設定
plt.xlabel('x')
plt.ylabel('密度')
plt.title(f'コーシー分布のシミュレーションと棄却域\nシミュレーションによる α = {alpha_simulated:.3f}')
plt.legend()
plt.grid(True)

# グラフを表示
plt.show()
シミュレーションによる第一種の過誤確率 α: 0.148

シミュレーション結果も手計算と一致しました。

2019 Q3(6)

唯一の不偏推定量であることを示しました。

 

コード

\hat{\theta} = \frac{n+1}{n} Yは、Yによるθの唯一の不偏推定量になります。またYはθの十分統計量ですから、完備十分統計量となります。よって、\hat{\theta} = \frac{n+1}{n} Yは、最小分散不偏推定量(UMVUE)となります。他のθの推定量2 \bar{X}と一緒にプロットし分散を確認します。

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
theta_true = 10  # 真のθ
n = 10  # サンプルサイズ
num_simulations = 10000  # シミュレーション回数

# UMVUE(最大値を使った推定量)の格納リスト
theta_hat_UMVUE_list = []

# 他の推定量(例えば、平均を使った推定量)の格納リスト
theta_hat_mean_list = []

# シミュレーション開始
for _ in range(num_simulations):
    # U(0, theta_true) から n 個の乱数を生成
    samples = np.random.uniform(0, theta_true, n)
    
    # 最大値 Y を使用して UMVUE を計算
    Y = np.max(samples)
    theta_hat_UMVUE = (n + 1) / n * Y
    theta_hat_UMVUE_list.append(theta_hat_UMVUE)
    
    # 平均を使用して他の推定量を計算
    theta_hat_mean = 2 * np.mean(samples)
    theta_hat_mean_list.append(theta_hat_mean)

# UMVUE 推定量の分散を計算
umvue_variance = np.var(theta_hat_UMVUE_list)

# 平均推定量の分散を計算
mean_variance = np.var(theta_hat_mean_list)

# 結果を表示
print(f"UMVUE(最大値を使った推定量)の分散: {umvue_variance}")
print(f"平均を使った推定量の分散: {mean_variance}")

# 推定量の分布をプロット
plt.hist(theta_hat_UMVUE_list, bins=30, alpha=0.7, label='UMVUE(最大値を使用)', color='blue', edgecolor='black', density=True)
plt.hist(theta_hat_mean_list, bins=30, alpha=0.7, label='平均を使用した推定量', color='red', edgecolor='black', density=True)
plt.title('UMVUE(最大値)と平均推定量の分布比較')
plt.xlabel('推定量の値')
plt.ylabel('密度')
plt.legend()
plt.grid(True)
plt.show()
UMVUE(最大値を使った推定量)の分散: 0.8726264971281747
平均を使った推定量の分散: 3.3430162264024963

\frac{n+1}{n} Yの分散は、2 \bar{X}の分散よりも小さいことが確認できました。

2019 Q3(5)

関数の期待値がゼロになる条件を求めました。

 

コード

u(Y)を与え、E[u(Y)]を求めるシミュレーションをします。

まずは適当なu(Y) = Y^2 - 2Yを与え、E[u(Y)]を求めてみます。

# 2019 Q3(5)  2024.9.21

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
theta_true = 10  # 真のθ
n = 15  # サンプルサイズ
num_simulations = 10000  # シミュレーション回数

# 仮の関数 u(Y) を定義(θと期待値に基づいた関数)
def u(Y, theta, n):
    # 例として、Y の関数として適当に定義(例えば u(Y) = Y^2 - 2Y)
    return Y**2 - 2*Y

# シミュレーションでの最大値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))

# 関数 u(Y) に基づく E[u(Y)] を計算
u_values = [u(y, theta_true, n) for y in max_values]
expected_u_Y = np.mean(u_values)

# 結果の表示
print(f"シミュレーションによる E[u(Y)]: {expected_u_Y}")

# u(Y) の値のヒストグラムを描画
plt.hist(u_values, bins=30, density=True, alpha=0.7, color='blue', edgecolor='black', label='u(Y) のシミュレーション結果')

# 期待値の線を追加
plt.axvline(expected_u_Y, color='red', linestyle='--', label=f'期待値 E[u(Y)] = {expected_u_Y:.4f}')

# グラフ設定
plt.title('関数 u(Y) の分布と期待値')
plt.xlabel('u(Y) の値')
plt.ylabel('密度')
plt.legend()
plt.grid(True)
plt.show()
シミュレーションによる E[u(Y)]: 69.66620056705366

E[u(Y)]=0にはなりませんでした。

次にu(Y) = 0を与え、E[u(Y)]を求めてみます。

# 2019 Q3(5)  2024.9.21

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
theta_true = 10  # 真のθ
n = 15  # サンプルサイズ
num_simulations = 10000  # シミュレーション回数

# 仮の関数 u(Y) を定義(θと期待値に基づいた関数)
def u(Y, theta, n):
    # 例: 0
    return 0

# シミュレーションでの最大値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))

# 関数 u(Y) に基づく E[u(Y)] を計算
u_values = [u(y, theta_true, n) for y in max_values]
expected_u_Y = np.mean(u_values)

# 結果の表示
print(f"シミュレーションによる E[u(Y)]: {expected_u_Y}")

# u(Y) の値のヒストグラムを描画
plt.hist(u_values, bins=30, density=True, alpha=0.7, color='blue', edgecolor='black', label='u(Y) のシミュレーション結果')

# 期待値の線を追加
plt.axvline(expected_u_Y, color='red', linestyle='--', label=f'期待値 E[u(Y)] = {expected_u_Y:.4f}')

# グラフ設定
plt.title('関数 u(Y) の分布と期待値')
plt.xlabel('u(Y) の値')
plt.ylabel('密度')
plt.legend()
plt.grid(True)
plt.show()
シミュレーションによる E[u(Y)]: 0.0

E[u(Y)]=0になりました。

次にu(Y) = Y - \frac{n}{n + 1} \thetaを与え、E[u(Y)]を求めてみます。

# 2019 Q3(5)  2024.9.21

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
theta_true = 10  # 真のθ
n = 15  # サンプルサイズ
num_simulations = 10000  # シミュレーション回数

# 仮の関数 u(Y) を定義(θと期待値に基づいた関数)
def u(Y, theta, n):
    # 例: Y - n/(n+1) * θ
    return Y - (n / (n + 1)) * theta

# シミュレーションでの最大値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))

# 関数 u(Y) に基づく E[u(Y)] を計算
u_values = [u(y, theta_true, n) for y in max_values]
expected_u_Y = np.mean(u_values)

# 結果の表示
print(f"シミュレーションによる E[u(Y)]: {expected_u_Y}")

# u(Y) の値のヒストグラムを描画
plt.hist(u_values, bins=30, density=True, alpha=0.7, color='blue', edgecolor='black', label='u(Y) のシミュレーション結果')

# 期待値の線を追加
plt.axvline(expected_u_Y, color='red', linestyle='--', label=f'期待値 E[u(Y)] = {expected_u_Y:.4f}')

# グラフ設定
plt.title('関数 u(Y) の分布と期待値')
plt.xlabel('u(Y) の値')
plt.ylabel('密度')
plt.legend()
plt.grid(True)
plt.show()
シミュレーションによる E[u(Y)]: -0.0069045066145089744

E[u(Y)]=0に近い値を取りました。u(Y) = 0でなくてもE[u(Y)]=0になるのはなぜでしょうか。おそらくこの問はu(Y)にθが含まれない前提になっているのでしょう。

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)とθの不偏推定量は理論値に近い値を取りました。

2019 Q3(3)

一様分布の最大値を条件とする条件付き同時確率密度関数の問題をやりました。

 

コード

一様分布の最大値を条件とする条件付き同時確率密度関数が

f(x_1, \dots, x_{n-1}, y | Y=y) = \frac{1}{y^{n-1}}

となるかをシミュレーションで確認しました。Yを1~10に変化させて、Y以下の乱数を発生させてカーネル密度推定で同時確率密度を計算します。

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

# パラメータ設定
n = 15  # サンプルサイズ
num_simulations = 10000  # シミュレーション回数
Y_values = np.arange(1, 11, 1)  # Y を 1 から 10 に変化させる

# シミュレーション結果を格納するリスト
simulated_f_values = []  # シミュレーションによる f の値
theoretical_f_values = []  # 理論値 f(x_1, ..., x_{n-1} | Y=y)

# シミュレーション開始
for Y in Y_values:
    f_sim_sum = 0  # f(x_1, ..., x_{n-1} | Y=y) の合計
    
    for _ in range(num_simulations):
        # Y より小さい範囲で n-1 個のサンプルを生成し、観測された最大値 Y を追加
        samples = np.random.uniform(0, Y, n-1)
        samples = np.append(samples, Y)  # 観測値としての Y を追加
        
        # カーネル密度推定を実行
        kde = gaussian_kde(samples)  # サンプルに基づいてKDEを実行
        
        # samples から Y (最後の要素) を除いた部分で密度を計算
        f_sim = np.prod(kde.evaluate(samples[:-1]))  # サンプルごとの確率密度を掛け合わせる
        
        # f(y) = 1 として扱う(条件付き分布)
        f_sim_sum += f_sim
    
    # シミュレーション結果の平均をリストに追加
    simulated_f_values.append(f_sim_sum / num_simulations)
    
    # 理論値 1 / Y^(n-1) を計算
    theoretical_f_values.append(1 / Y**(n-1))

# グラフ描画
plt.plot(Y_values, simulated_f_values, 'bo-', label='シミュレーション結果')
plt.plot(Y_values, theoretical_f_values, 'r--', label='理論値 1/Y^(n-1)')

# 縦軸を対数スケールに変更
plt.yscale('log')

# グラフ設定
plt.title(f'f(x1, ..., xn-1 | Y=y) のシミュレーション結果と理論値 (n={n})')
plt.xlabel('Y の値')
plt.ylabel('f(x1, ..., xn-1 | Y=y)')
plt.legend()
plt.grid(True)
plt.show()

概ね理論値と一致しました。カーネル密度推定の誤差により、シミュレーションの値は理論値より少し大きくなっていると思われます。

2019 Q3(2)

一様分布の最大値をとる確率変数の確率密度関数を求めました。

 

コード

数式を使って確率密度関数g(y)を求めます。

# 2019 Q3(2)  2024.9.18

import sympy as sp

# 変数を定義
y, theta, n = sp.symbols('y theta n')

# 累積分布関数 F_Y(y) を定義
F_Y = (y / theta)**n  # F_Y(y) = (y / θ)^n

# 確率密度関数 g(y) を F_Y(y) の y に関する微分として計算
g_y = sp.diff(F_Y, y)

# 簡略化
g_y_simplified = sp.simplify(g_y)

# 結果を表示
display(g_y_simplified)

形は少し違いますが、正しいです。

得られた確率密度関数g(y)と数値シミュレーションによる分布と比較します。

import numpy as np
import matplotlib.pyplot as plt

# パラメータ
n = 10  # サンプルサイズ(n個の一様分布からのサンプル)
theta = 10  # 一様分布の上限
num_simulations = 10000  # シミュレーション回数

# 最大値 Y のリスト
Y_values = []

# シミュレーション
for _ in range(num_simulations):
    samples = np.random.uniform(0, theta, n)  # (0, θ)の範囲からn個のサンプルを生成
    Y = np.max(samples)  # 最大値を計算
    Y_values.append(Y)

# ヒストグラムを描画
plt.hist(Y_values, bins=30, density=True, alpha=0.7, color='blue', edgecolor='black', label='シミュレーション結果')

# 理論上の確率密度関数を描画
y = np.linspace(0, theta, 1000)
g_y = (n / theta**n) * y**(n-1)
plt.plot(y, g_y, color='red', label='理論的な分布')

# グラフの設定
plt.title(f'最大値 Y の分布と理論的分布 (n={n}, θ={theta})')
plt.xlabel('最大値 Y')
plt.ylabel('確率密度')
plt.legend()
plt.show()

確率密度関数g(y)と数値シミュレーションによる分布は一致することが確認できました。