2019 Q2(4)
損失関数の期待値を導出しそれが最小となるパラメータαを求める問題をやりました。
コード
数値シミュレーションによりR(α)が最小値になるαを見つけます。
# 2019 Q2(4) 2024.9.16
import numpy as np
import matplotlib.pyplot as plt
# 損失関数 R(α) の定義
def R(alpha):
return alpha + 2 / alpha - 2
# α の範囲を設定 (0.1 から 3 まで、細かいステップで計算)
alpha_range = np.linspace(0.1, 3, 1000)
# 損失関数 R(α) の計算
R_values = R(alpha_range)
# 最小値を取る α の確認 (numpy の argmin を使って最小値を見つける)
min_index = np.argmin(R_values)
min_alpha = alpha_range[min_index]
min_R = R_values[min_index]
# グラフのプロット
plt.figure(figsize=(8, 6))
plt.plot(alpha_range, R_values, label='R(α)')
plt.axvline(np.sqrt(2), color='red', linestyle='--', label=r'$\alpha = \sqrt{2}$')
plt.scatter(min_alpha, min_R, color='red', zorder=5, label=f'最小値: α={min_alpha:.2f}, R(α)={min_R:.2f}')
plt.xlabel(r'$\alpha$')
plt.ylabel(r'$R(\alpha)$')
plt.title(r'$R(\alpha)$ の最小値確認 (数値シミュレーション)')
plt.legend()
plt.grid(True)
plt.show()
# 最小値の表示
print(f'数値シミュレーションでの最小値 α: {min_alpha}, R(α): {min_R}')
数値シミュレーションでの最小値 α: 1.4150150150150151, R(α): 0.8284275786822475
手計算と近い値になりました。
極値付近が少し平らなので、αの範囲を1.0~2.0に変更し拡大してみます。
import numpy as np
import matplotlib.pyplot as plt
# 損失関数 R(α) の定義
def R(alpha):
return alpha + 2 / alpha - 2
# α の範囲を 1 から 2 に設定
alpha_range = np.linspace(1, 2, 1000)
# 損失関数 R(α) の計算
R_values = R(alpha_range)
# 最小値を取る α の確認 (numpy の argmin を使って最小値を見つける)
min_index = np.argmin(R_values)
min_alpha = alpha_range[min_index]
min_R = R_values[min_index]
# グラフのプロット
plt.figure(figsize=(8, 6))
plt.plot(alpha_range, R_values, label='R(α)')
plt.axvline(np.sqrt(2), color='red', linestyle='--', label=r'$\alpha = \sqrt{2}$')
plt.scatter(min_alpha, min_R, color='red', zorder=5, label=f'最小値: α={min_alpha:.4f}, R(α)={min_R:.4f}')
plt.xlabel(r'$\alpha$')
plt.ylabel(r'$R(\alpha)$')
plt.title(r'$R(\alpha)$ の最小値確認 (αの範囲: 1 から 2)')
plt.legend()
plt.grid(True)
plt.show()
# 最小値の表示
print(f'数値シミュレーションでの最小値 α: {min_alpha}, R(α): {min_R}')
数値シミュレーションでの最小値 α: 1.4144144144144144, R(α): 0.8284271532679175
最小値が取れていることが分かりました。
2019 Q2(3)
(2)で求めた確率密度関数において確率変数Uの逆数の期待値を求めました。
コード
数式を使って1/Uの期待値を求めます。
# 2019 Q2(3) 2024.9.15
import sympy as sp
# 変数の定義
u, lambda_value = sp.symbols('u lambda', positive=True, real=True)
# ガンマ分布の確率密度関数 g(u) = λ^2 * u * exp(-λ * u)
g_u = lambda_value**2 * u * sp.exp(-lambda_value * u)
# 1/U の期待値の積分
integrand = (1/u) * g_u # 1/u * g(u)
# 積分の実行
expectation = sp.integrate(integrand, (u, 0, sp.oo))
# 結果の簡略化
expectation_simplified = sp.simplify(expectation)
# 結果の表示
display(expectation_simplified)
手計算と一致します。
次に、数値シミュレーションを行います。
import numpy as np
import matplotlib.pyplot as plt
# パラメータ設定
lambda_value = 2 # λ=2
num_samples = 100000 # サンプルの数
# 指数分布に従う乱数 X1, X2 を生成
X1 = np.random.exponential(scale=1/lambda_value, size=num_samples)
X2 = np.random.exponential(scale=1/lambda_value, size=num_samples)
# U = X1 + X2 の計算
U = X1 + X2
# 1/U の計算
inv_U = 1 / U
# シミュレーションの結果からの平均値
simulation_result = np.mean(inv_U)
# 理論値 E[1/U] = λ
theoretical_value = lambda_value
# 結果の表示
print(f'シミュレーションによる E[1/U]: {simulation_result}')
print(f'理論値 E[1/U]: {theoretical_value}')
# ヒストグラムの描画
plt.hist(inv_U, bins=50, density=True, alpha=0.7, label='シミュレーション結果')
plt.axvline(theoretical_value, color='r', linestyle='--', label=f'理論値 E[1/U] = {theoretical_value}')
plt.xlabel('1/U')
plt.ylabel('確率密度')
plt.title('シミュレーションによる 1/U の分布と理論値の比較')
plt.legend()
plt.show()
シミュレーションによる E[1/U]: 1.9928303711222939
理論値 E[1/U]: 2
シミュレーションによりE[1/U]は理論値に近づきました。分散は発散するためヒストグラムにすると大きなバラつきが見られます。
2019 Q2(2)
独立な指数分布に従う2変数の和の確率密度関数を求めました。
コード
数式を使ってUの確率密度関数を求めます。
# 2019 Q2(2) 2024.9.14
import sympy as sp
# 変数の定義
x, u, lambda_value = sp.symbols('x u lambda', positive=True, real=True)
# 指数分布の確率密度関数 (PDF) f(x)
f_x = lambda_value * sp.exp(-lambda_value * x)
# 畳み込み積分を計算して、g(u) = f(x) * f(u - x) の形式にする
g_u = sp.integrate(f_x * lambda_value * sp.exp(-lambda_value * (u - x)), (x, 0, u))
# 簡略化
g_u_simplified = sp.simplify(g_u)
# 結果の表示
display(g_u_simplified)
X1 (X2)と、Uの確率密度関数を重ねて描画してみます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma
# パラメータ設定
lambda_value = 2 # λ=2
x_range = np.linspace(0, 5, 1000) # X1, X2 の範囲
# X1, X2 の指数分布の確率密度関数 (PDF) - 同じなので1つにまとめる
pdf_X1_X2 = lambda_value * np.exp(-lambda_value * x_range)
# U = X1 + X2 は形状パラメータk=2、スケールパラメータθ=1/λのガンマ分布
u_range = np.linspace(0, 10, 1000)
pdf_U = gamma.pdf(u_range, a=2, scale=1/lambda_value)
# グラフのプロット
plt.figure(figsize=(8, 6))
# X1 (X2) の指数分布の描画
plt.plot(x_range, pdf_X1_X2, label='X1 (X2) (指数分布)', color='blue', linestyle='--')
# U のガンマ分布の描画
plt.plot(u_range, pdf_U, label='U = X1 + X2 (ガンマ分布)', color='red')
# グラフの装飾
plt.title('X1 (X2) と U の確率密度関数の比較')
plt.xlabel('値')
plt.ylabel('確率密度')
plt.legend()
# グラフの表示
plt.show()
グラフの形状から、U = X1 + X2の関係は想像しづらいですね。
ガンマ分布の形状パラメータを1~2に変化させて、X1 (X2)からUへの変化を確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma
# パラメータ設定
lambda_value = 2 # λ=2
x_range = np.linspace(0, 5, 1000) # X1, X2, U の範囲を0~5に設定
# ガンマ分布の形状パラメータを1から2まで0.1ずつ変化させる
shape_params = np.arange(1, 2.1, 0.1) # 1から2までの形状パラメータを0.1ステップで変化
# グラフのプロット
plt.figure(figsize=(8, 6))
# 形状パラメータが1から2に変化するガンマ分布の描画
for k in shape_params:
gamma_pdf = gamma.pdf(x_range, a=k, scale=1/lambda_value)
plt.plot(x_range, gamma_pdf, label=f'形状パラメータ k={k:.1f}')
# X1 (X2) の指数分布の描画 - 同じなので1つにまとめる
pdf_X1_X2 = lambda_value * np.exp(-lambda_value * x_range)
plt.plot(x_range, pdf_X1_X2, label='X1 (X2) (指数分布)', color='blue', linestyle='--')
# グラフの装飾
plt.title('指数分布からガンマ分布への変化 (形状パラメータの0.1ステップ変化)')
plt.xlabel('値')
plt.ylabel('確率密度')
plt.legend()
# グラフの表示
plt.show()
X1 (X2)からUへの変化が可視化できました。
次に、X1 (X2)と、Uの累積分布関数を重ねて描画してみます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import expon, gamma
# パラメータ設定
lambda_value = 2 # λ=2
x_range = np.linspace(0, 5, 1000) # X1, X2 の範囲
# X1, X2 の指数分布の累積分布関数 (CDF) - 同じなので1つにまとめる
cdf_X1_X2 = expon.cdf(x_range, scale=1/lambda_value)
# U = X1 + X2 は形状パラメータk=2、スケールパラメータθ=1/λのガンマ分布
u_range = np.linspace(0, 10, 1000)
cdf_U = gamma.cdf(u_range, a=2, scale=1/lambda_value)
# グラフのプロット
plt.figure(figsize=(8, 6))
# X1 (X2) の指数分布のCDFの描画
plt.plot(x_range, cdf_X1_X2, label='X1 (X2) (指数分布 CDF)', color='blue', linestyle='--')
# U のガンマ分布のCDFの描画
plt.plot(u_range, cdf_U, label='U = X1 + X2 (ガンマ分布 CDF)', color='red')
# グラフの装飾
plt.title('X1 (X2) と U の累積分布関数 (CDF) の比較')
plt.xlabel('値')
plt.ylabel('累積確率')
plt.legend()
# グラフの表示
plt.show()
グラフの形状から、U = X1 + X2の関係を想像しやすくなりました。
2019 Q2(1)
独立な指数分布に従う2変数の和の期待値を求めました。
コード
数式を使ってUの期待値を求めます
import sympy as sp
# 変数の定義
x = sp.Symbol('x', positive=True)
lambda_value = sp.Symbol('lambda', positive=True)
# 指数分布の確率密度関数 (PDF)
f_x = lambda_value * sp.exp(-lambda_value * x)
# 期待値の式を計算 (積分)
E_X1 = sp.integrate(x * f_x, (x, 0, sp.oo))
# 期待値
E_X1_simplified = sp.simplify(E_X1)
# X1 + X2 の場合 (独立性より)
E_U = 2 * E_X1_simplified
# 結果の表示
display(E_U)
手計算と一致します
次にUに従う乱数を生成し分布を確認します。また形状パラメータ2のガンマ分布と重ねて描画します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma
# パラメータ設定
lambda_value = 2 # 例としてλ=2
num_samples = 1000 # サンプルの数
# 指数分布に従う乱数の生成
X1 = np.random.exponential(scale=1/lambda_value, size=num_samples)
X2 = np.random.exponential(scale=1/lambda_value, size=num_samples)
# U = X1 + X2 の計算
U = X1 + X2
# 理論的な期待値
expected_value = 2 / lambda_value
# シミュレーション結果のヒストグラムのプロット (確率密度にスケール)
plt.hist(U, bins=30, density=True, alpha=0.7, label='シミュレーション結果')
# ガンマ分布 (形状パラメータ k=2, スケールパラメータ θ=1/λ) の確率密度関数 (PDF) を計算
k = 2 # 形状パラメータ (k = 2)
theta = 1 / lambda_value # スケールパラメータ (θ = 1/λ)
# ガンマ分布に基づくx軸の範囲
x = np.linspace(0, max(U), 1000)
# ガンマ分布の確率密度関数 (PDF) を計算
gamma_pdf = gamma.pdf(x, a=k, scale=theta)
# ガンマ分布の折れ線グラフをプロット
plt.plot(x, gamma_pdf, 'r-', label=f'ガンマ分布 (k={k}, θ={theta:.2f})')
# グラフの装飾
plt.axvline(expected_value, color='g', linestyle='--', label=f'理論値 = {expected_value:.2f}')
plt.xlabel('U = X1 + X2')
plt.ylabel('確率密度')
plt.title('Uの分布とガンマ分布の比較')
plt.legend()
plt.show()
# 実際のサンプル平均値と理論的期待値の比較
print(f'理論的な期待値: {expected_value}')
print(f'シミュレーションからの平均値: {np.mean(U)}')
理論的な期待値: 1.0
シミュレーションからの平均値: 0.9902440726761544
Uの分布と形状パラメータ2のガンマ分布は重なりました。Uはガンマ分布に従うことが分かります。
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
2024 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
2022 Q4(4)
対数で変数変換した場合の分布を求めました。
コード
数式を使った計算
# 2022 Q4(4) 2024.8.11
import sympy as sp
# シンボリック変数の定義
t, gamma = sp.symbols('t gamma', positive=True)
x = sp.symbols('x', real=True, positive=True)
# 累積分布関数 F_X(x) の定義
F_X = sp.Piecewise((0, x <= 1), (1 - x**(-1/gamma), x > 1))
# F_T(t) の計算
F_T = F_X.subs(x, sp.exp(gamma * t))
# f_T(t) の計算(F_T(t) を t で微分)
f_T = sp.diff(F_T, t)
# 結果の表示
display(F_T, f_T)
# LaTeXで表示できなときは、こちら
#sp.pprint(F_T, use_unicode=True)
#sp.pprint(f_T, use_unicode=True)
シミュレーションによる計算
# 2022 Q4(4) 2024.8.11
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
# パラメータ設定
gamma = 1.0 # γの値(例えば1)
num_samples = 10000 # サンプルサイズ
# 累積分布関数 F(x) に従う乱数を生成
X = (np.random.uniform(size=num_samples))**(-gamma)
# T = (1/γ) * log(X) の計算
T = (1/gamma) * np.log(X)
# ヒストグラムをプロット
plt.hist(T, bins=50, density=True, alpha=0.6, color='g', label="シミュレーションによる T")
# 理論的な指数分布の密度関数をプロット
x = np.linspace(0, 8, 100)
pdf = stats.expon.pdf(x) # 指数分布(λ=1)の密度関数
plt.plot(x, pdf, 'r-', lw=2, label='指数分布の理論値')
# グラフの設定(日本語でキャプションとラベルを設定)
plt.xlabel('Tの値')
plt.ylabel('確率密度')
plt.title(r'$T = \frac{1}{\gamma} \log X$ のヒストグラムと指数分布の比較')
plt.legend()
plt.show()