ホーム » 分布 (ページ 7)

「分布」カテゴリーアーカイブ

2019 Q5(3)

投稿日:2024年3月11日　最終更新日:2024年9月29日

前問の事後確率密度関数の最大値を取る期待値μを求めました。

コード

対数事後分布 h(μ)と、その第1項（観測データに基づく項）および第2項（事前分布に基づく項）のグラフを描画します。パラメータ ξ を 1 または 5、スケールパラメータ λ を 25 または 50 として、それぞれの組み合わせで描画します。

# 2019 Q5(3)  2024.9.29

import numpy as np
import matplotlib.pyplot as plt

# 対数事後分布の第1項: 観測データに基づく項
def h_term1(mu, y):
    n = len(y)
    y_mean = np.mean(y)
    return -0.5 * n * (mu - y_mean)**2

# 対数事後分布の第2項: 事前分布に基づく項
def h_term2(mu, lambd, xi):
    return -lambd * np.abs(mu - xi)

# 対数事後分布 h(μ) の定義
def h(mu, y, lambd, xi):
    return h_term1(mu, y) + h_term2(mu, lambd, xi)

# μハット（事後分布のモード）を計算
def mu_hat(y_mean, lambd, xi, n):
    if y_mean > xi:
        return max(y_mean - lambd / n, xi)
    elif y_mean < xi:
        return min(y_mean + lambd / n, xi)
    else:
        return xi

# パラメータ設定のリスト (xi, lambd) の組み合わせ
params = [(1, 25), (1, 50), (5, 25), (5, 50)]
n_samples = 20  # サンプル数
mu_values = np.linspace(0, 6, 500)  # μの範囲

# グラフを4つ表示するための設定
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

for i, (xi, lambd) in enumerate(params):
    # 観測データを生成して、平均を3に調整
    np.random.seed(43)
    y = np.random.normal(loc=3, scale=1, size=n_samples)
    y_mean = np.mean(y)
    
    # μハットを計算
    mu_hat_value = mu_hat(y_mean, lambd, xi, n_samples)

    # 各項の計算
    h_values = h(mu_values, y, lambd, xi)
    h_term1_values = h_term1(mu_values, y)
    h_term2_values = h_term2(mu_values, lambd, xi)

    # サブプロットに描画
    ax = axes[i//2, i%2]
    ax.plot(mu_values, h_values, label='h(μ)', color='orange')
    ax.plot(mu_values, h_term1_values, label='第1項: -n/2 (μ - ȳ)²', color='blue')
    ax.plot(mu_values, h_term2_values, label='第2項: -λ|μ - ξ|', color='green')
    
    # μハットを縦線でプロット
    ax.axvline(mu_hat_value, color='red', linestyle='--', label=f'μハット = {mu_hat_value:.2f}')

    # グラフのタイトル
    ax.set_title(f'ξ = {xi}, λ = {lambd}, ȳ = {y_mean:.2f}, n = {n_samples}')
    ax.set_xlabel('μ')
    ax.set_ylabel('値')
    ax.legend()
    ax.grid(True)

# グラフを描画
plt.tight_layout()
plt.show()

対数事後分布 h(μ)が最大値を取るμハットが

であることが確認できました。

2019 Q5(2)

投稿日:2024年3月10日　最終更新日:2024年9月28日

ラプラス分布に従う確率変数μを母平均とする正規分布の観測値からμの事後確率密度関数を求めました。

コード

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

# 2019 Q5(2)  2024.9.28

import sympy as sp

# シンボリック変数の定義
mu, xi, lambd, n = sp.symbols('mu xi lambda n')
y = sp.IndexedBase('y')
i = sp.symbols('i', integer=True)

# ラプラス分布の事前分布
laplace_prior = (lambd / 2) * sp.exp(-lambd * sp.Abs(mu - xi))

# 正規分布の尤度関数
likelihood = (1 / (2 * sp.pi)**(n / 2)) * sp.exp(-sp.Sum((y[i] - mu)**2 / 2, (i, 1, n)))

# 事後分布 (比例定数は無視)
posterior = laplace_prior * likelihood

# 結果を簡略化
posterior_simplified = sp.simplify(posterior)

# 結果を表示
posterior_simplified

形は少し違いますが手計算の結果と一致しました。

次に、事前分布と事後分布の確率密度関数をプロットします。比較しやすいように事後分布は正規化します。また事前分布がない場合として尤度関数も重ねて描画してみましょう。

# 2019 Q5(2)  2024.9.28

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simps

# 事前分布（ラプラス分布）の定義
def laplace_prior(mu, xi, lambd):
    return (lambd / 2) * np.exp(-lambd * np.abs(mu - xi))

# 事後分布（正規化定数は省略）
def posterior_example(mu, y, xi, lambd):
    n = len(y)
    y_mean = np.mean(y)
    sum_y = np.sum((y - y_mean)**2)
    
    # 解答例の事後分布
    return (lambd / (2 * (2 * np.pi)**(n / 2))) * np.exp(-0.5 * (sum_y + n * (mu - y_mean)**2) - lambd * np.abs(mu - xi))

# 尤度関数（事前分布がない場合）
def likelihood(mu, y):
    n = len(y)
    y_mean = np.mean(y)
    return np.exp(-0.5 * n * (mu - y_mean)**2)

# パラメータ設定
xi = 0  # ラプラス分布の中央値
lambd = 1  # スケールパラメータ
n_samples = 20  # サンプル数

# 観測データを生成
np.random.seed(43)
y = np.random.normal(loc=0, scale=1, size=n_samples)

# μの範囲を設定
mu_values = np.linspace(-3, 3, 500)

# 事前分布の計算
prior_values = laplace_prior(mu_values, xi, lambd)

# 事後分布の計算（事前分布あり）
posterior_values_example = posterior_example(mu_values, y, xi, lambd)

# 尤度関数の計算（事前分布なし）
likelihood_values = likelihood(mu_values, y)

# 事後分布と尤度関数を正規化
area_under_curve_example = simps(posterior_values_example, mu_values)  # 数値積分
posterior_values_example /= area_under_curve_example  # 正規化

area_under_likelihood = simps(likelihood_values, mu_values)  # 数値積分
likelihood_values /= area_under_likelihood  # 正規化

# グラフをプロット
plt.plot(mu_values, prior_values, label='事前分布 (ラプラス分布)', color='blue')
plt.plot(mu_values, posterior_values_example, label='事後分布', color='purple')
plt.plot(mu_values, likelihood_values, label='尤度関数（事前分布なし）', color='green', linestyle='--')
plt.title('事前分布、事後分布、および事前分布がない場合の尤度関数の比較')
plt.xlabel('μ')
plt.ylabel('確率密度')
plt.legend()
plt.show()

事後分布は事前分布と尤度関数の中間的な形をとっていることが分かります。

2019 Q5(1)-2

投稿日:2024年3月9日　最終更新日:2024年9月28日

ラプラス分布の分散を求めました。

コード

ξ=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 Q5(1)-1

投稿日:2024年3月8日　最終更新日:2024年9月28日

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

コード

ξ=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)

投稿日:2024年3月7日　最終更新日:2024年9月26日

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

コード

棄却域を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)

投稿日:2024年3月6日　最終更新日:2024年9月25日

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

コード

尤度比λ(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)

投稿日:2024年3月5日　最終更新日:2024年9月24日

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

コード

数式を使って検出力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)

投稿日:2024年3月4日　最終更新日:2024年9月23日

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

コード

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

# 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(5)

投稿日:2024年3月2日　最終更新日:2024年9月21日

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

コード

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)

投稿日:2024年3月1日　最終更新日:2025年1月31日

一様分布に従う複数の確率変数がとる最大値の期待値と、その最大値から上限θの不偏推定量を求めました。

コード

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

「分布」カテゴリーアーカイブ

カテゴリー

投稿一覧