2018 Q5(1)
一様分布の順序統計量の確率密度関数と期待値を求めました。
コード
シミュレーションにより順序統計量Y1とY3の分布と期待値が理論値と一致するか確認します。
# 2018 Q5(1) 2024.10.21
import numpy as np
import matplotlib.pyplot as plt
# 定義域
y_values = np.linspace(0, 1, 500)
# 確率密度関数の定義
f1_y = 3 * (1 - y_values) ** 2
f3_y = 3 * y_values ** 2
# シミュレーションの設定
n_samples = 10000 # サンプル数
# 一様分布から3つの変数をサンプリング
samples = np.random.uniform(0, 1, (n_samples, 3))
# 順序統計量を求める
Y1_samples = np.min(samples, axis=1) # 最小値
Y3_samples = np.max(samples, axis=1) # 最大値
# 期待値の計算
expected_Y1 = np.mean(Y1_samples)
expected_Y3 = np.mean(Y3_samples)
# 理論値との比較
theoretical_Y1 = 1/4
theoretical_Y3 = 3/4
# 期待値を表示
print(f"シミュレーションによる E[Y1](最小値の期待値): {expected_Y1:.4f}, 理論値: {theoretical_Y1:.4f}")
print(f"シミュレーションによる E[Y3](最大値の期待値): {expected_Y3:.4f}, 理論値: {theoretical_Y3:.4f}")
# ヒストグラムの作成
plt.figure(figsize=(8, 6))
# 理論密度関数のプロット
plt.plot(y_values, f1_y, label='$f_1(y) = 3(1-y)^2$', color='blue')
plt.plot(y_values, f3_y, label='$f_3(y) = 3y^2$', color='green')
# サンプルのヒストグラムを追加 (確率密度として正規化)
plt.hist(Y1_samples, bins=50, density=True, alpha=0.5, color='blue', label='Y1 samples (min)')
plt.hist(Y3_samples, bins=50, density=True, alpha=0.5, color='green', label='Y3 samples (max)')
# グラフの設定
plt.title('確率密度関数 $f_1(y)$ と $f_3(y)$ (シミュレーションと理論)', fontsize=14)
plt.xlabel('y', fontsize=12)
plt.ylabel('密度', fontsize=12)
plt.legend(loc='upper right', fontsize=10)
plt.grid(True)
# グラフの表示
plt.show()
シミュレーションによる E[Y1](最小値の期待値): 0.2499, 理論値: 0.2500
シミュレーションによる E[Y3](最大値の期待値): 0.7487, 理論値: 0.7500
シミュレーションによる順序統計量Y1とY3の分布と期待値が理論値と一致しました。
2019 Q3(5)
関数の期待値がゼロになる条件を求めました。
コード
u(Y)を与え、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
にはなりませんでした。
次にを与え、を求めてみます。
# 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
になりました。
次にを与え、を求めてみます。
# 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
に近い値を取りました。でなくてもになるのはなぜでしょうか。おそらくこの問はにθが含まれない前提になっているのでしょう。
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)
一様分布の最大値を条件とする条件付き同時確率密度関数の問題をやりました。
コード
一様分布の最大値を条件とする条件付き同時確率密度関数が
となるかをシミュレーションで確認しました。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)と数値シミュレーションによる分布は一致することが確認できました。
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はθに関する十分統計量であることが確認できました。
2021 Q1(4)
独立ではない2つの確率変数の積の期待値を求めました。
コード
数式を使った計算
# 2021 Q1(4) 2024.8.20
import sympy as sp
# シンボリック変数の定義
x, y = sp.symbols('x y')
# h(x) の導出
h_x = sp.integrate(sp.exp(-x), (x, 0, x))
# XY の期待値の計算
# h(x) = 1 - exp(-x) を使って E[X * h(X)] を計算する
integral_part1 = sp.integrate(x * h_x * sp.exp(-x), (x, 0, sp.oo))
E_XY_simplified = sp.simplify(integral_part1)
# 結果の表示
display(sp.Eq(sp.Symbol('h(x)'), h_x))
display(sp.Eq(sp.Symbol('E[XY]'), E_XY_simplified))
指数分布に従う乱数Xを生成し、Y=h(X)Yが一様分布に従うか確認
import numpy as np
import matplotlib.pyplot as plt
# シミュレーションの設定
num_samples = 100000 # サンプル数
# 指数分布に従う乱数 X を生成
X_samples = np.random.exponential(scale=1.0, size=num_samples)
# 関数 Y = h(X) = 1 - exp(-X) を適用して Y を計算
Y_samples = 1 - np.exp(-X_samples)
# ヒストグラムのプロット
plt.figure(figsize=(8, 5))
plt.hist(Y_samples, bins=50, density=True, alpha=0.6, color='g', label='シミュレーションされた Y=h(X)')
plt.xlabel('y')
plt.ylabel('密度')
plt.title('Y = h(X) の分布 (シミュレーション結果)')
plt.legend()
plt.grid(True)
plt.show()
シミュレーションによるE[XY]の計算
import numpy as np
# シミュレーションの設定
num_samples = 100000 # サンプル数
# 指数分布に従う乱数 X を生成
X_samples = np.random.exponential(scale=1.0, size=num_samples)
# 関数 Y = h(X) = 1 - exp(-X) を適用して Y を計算
Y_samples = 1 - np.exp(-X_samples)
# XY の期待値を計算
XY_samples = X_samples * Y_samples
E_XY_simulation = np.mean(XY_samples)
# 結果の表示
print(f"シミュレーションによる E[XY] = {E_XY_simulation}")
print(f"理論値 E[XY] = 3/4 = {3/4}")
シミュレーションによる E[XY] = 0.7473462073557277
理論値 E[XY] = 3/4 = 0.75
2021 Q1(3)
変数変換を用いて二つの分布の和の確率密度関数を求めました。
コード
数式を使った計算
# 2024 Q1(3) 2024.8.19
import sympy as sp
# シンボリック変数の定義
x, z = sp.symbols('x z')
# 確率密度関数の定義
f_X = sp.exp(-x) # f_X(x) = e^(-x)
f_Y = sp.Piecewise((1, (z - x >= 0) & (z - x <= 1)), (0, True)) # f_Y(z - x) for 0 < z - x < 1
# 一般的な積分の設定
f_Z_integral = sp.integrate(f_X * f_Y, (x, 0, z))
# 結果の簡略化
f_Z_general = sp.simplify(f_Z_integral)
# 結果の表示
display(f_Z_general)
簡単のため、zの範囲を指定して再度計算
# 2024 Q1(3) 2024.8.19
import sympy as sp
# シンボリック変数の定義
x, z = sp.symbols('x z')
# 確率密度関数の定義
f_X = sp.exp(-x) # f_X(x) = e^(-x)
f_Y = 1 # f_Y(y) = 1 (0 < y < 1)
# 各範囲での計算
# 1. z <= 0 の場合
f_Z_1 = 0 # z <= 0 の場合は f_Z(z) = 0
# 2. 0 < z <= 1 の場合
f_Z_2_integral = sp.integrate(f_X * f_Y, (x, 0, z))
f_Z_2 = sp.simplify(f_Z_2_integral)
# 3. z > 1 の場合
f_Z_3_integral = sp.integrate(f_X * f_Y, (x, z-1, z))
f_Z_3 = sp.simplify(f_Z_3_integral)
# 結果の表示
print(f"f_Z(z) for z <= 0:")
display(f_Z_1)
print(f"f_Z(z) for 0 < z <= 1:")
display(f_Z_2)
print(f"f_Z(z) for z > 1:")
display(f_Z_3)
プロット
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
# 確率密度関数 f_X(x) と f_Y(y)
def f_X(x):
return np.exp(-x) if x > 0 else 0
def f_Y(y):
return 1 if 0 < y < 1 else 0
# Z = X + Y の確率密度関数 f_Z(z) の計算
def f_Z(z):
integrand = lambda x: f_X(x) * f_Y(z - x)
return quad(integrand, max(0, z-1), z)[0] #積分範囲を限定する場合
#return quad(integrand, 0, z, epsabs=1e-8, epsrel=1e-8)[0] #積分範囲を広くする場合
# z の範囲を設定し、f_Z(z) を計算
z_values = np.linspace(0, 5, 1000) # 増加させたサンプル数
f_Z_values = np.array([f_Z(z) for z in z_values])
# グラフの描画
plt.plot(z_values, f_Z_values, label='f_Z(z) = X + Y')
plt.xlabel('z')
plt.ylabel('確率密度')
plt.title('X + Y の確率密度関数')
plt.legend()
plt.grid(True)
plt.show()
2021 Q1(1)(2)
指数分布と一様分布の期待値と、それらが独立な場合の積の期待値を求めました。
コード
数式を使った計算
# 2021 Q1(1)(2) 2024.8.18
import numpy as np
from scipy.integrate import quad
# Xの期待値の計算
def fx(x):
return x * np.exp(-x)
E_X, _ = quad(fx, 0, np.inf)
# Yの期待値の計算
def fy(y):
return y
E_Y, _ = quad(fy, 0, 1)
# XYの期待値の計算(独立の場合)
E_XY = E_X * E_Y
print(f"E[X] = {E_X}")
print(f"E[Y] = {E_Y}")
print(f"E[XY] (独立の場合) = {E_XY}")
E[X] = 0.9999999999999998
E[Y] = 0.5
E[XY] (独立の場合) = 0.4999999999999999
シミュレーションによる計算
# 2021 Q1(1)(2) 2024.8.18
import numpy as np
# シミュレーションの設定
num_samples = 1000000 # サンプル数
# Xの乱数生成(指数分布)
X_samples = np.random.exponential(scale=1.0, size=num_samples)
# Yの乱数生成(一様分布)
Y_samples = np.random.uniform(0, 1, size=num_samples)
# (1) XとYの期待値の計算
E_X_simulation = np.mean(X_samples)
E_Y_simulation = np.mean(Y_samples)
# (2) XYの期待値の計算(独立な場合)
E_XY_simulation = np.mean(X_samples * Y_samples)
print(f"E[X] = {E_X_simulation}")
print(f"E[Y] = {E_Y_simulation}")
print(f"E[XY] = {E_XY_simulation}")
E[X] = 1.0001064255963075
E[Y] = 0.4998284579940602
E[XY] = 0.4998038319582688