青空統計学

# 2022 Q4(3)  2024.8.10

import sympy as sp

# シンボリック変数の定義
gamma = sp.Symbol('gamma', positive=True)
n = 10  # 具体的なサンプル数を設定
X = sp.symbols('X1:%d' % (n+1), positive=True)  # X1, X2, ..., Xn を定義

# 尤度関数の定義
L_gamma = sp.prod((1/gamma) * X[i]**(-(1/gamma + 1)) for i in range(n))

# 対数尤度関数の定義
log_L_gamma = sp.log(L_gamma)

# 對数尤度関数を γ で微分
diff_log_L_gamma = sp.diff(log_L_gamma, gamma)

# 微分した結果を 0 に等しいとおいて解く
gamma_hat_eq = sp.Eq(diff_log_L_gamma, 0)
gamma_hat = sp.solve(gamma_hat_eq, gamma)

# 結果の表示
sp.pprint(gamma_hat[0], use_unicode=True)

# 2022 Q4(3)  2024.8.10

import numpy as np

# 観測データ（例としてランダムに生成）
n = 1000  # データ点の数
gamma_true = 0.5  # 実際の γ (例として)

# CDFの逆関数
X = (np.random.uniform(size=n))**(-gamma_true)

# 最尤推定量 gamma_hat} の計算
log_X = np.log(X)
gamma_hat = np.mean(log_X)

print(f"観測(ランダム)データからの最尤推定量 gamma_hat: {gamma_hat}")

観測(ランダム)データからの最尤推定量 gamma_hat: 0.4984974482026361

# 2022 Q4(2)  2024.8.9

from sympy import init_printing
import sympy as sp

# これでTexの表示ができる
init_printing()

# シンボリック変数の定義
x, gamma = sp.symbols('x gamma', positive=True)

# 確率密度関数 f(x) の定義
f_x = (1/gamma) * x**(-(1/gamma + 1))

# 期待値 E[X] の計算
E_X = sp.integrate(x * f_x, (x, 1, sp.oo))
E_X_simplified = sp.simplify(E_X)

# E[X^2] の計算
E_X2 = sp.integrate(x**2 * f_x, (x, 1, sp.oo))
E_X2_simplified = sp.simplify(E_X2)

# 分散 V[X] の計算
V_X = E_X2_simplified - E_X_simplified**2

# 分散をさらに簡略化
V_X_simplified = sp.simplify(V_X)
V_X_factored = sp.factor(V_X_simplified)

# 結果を表示
E_X_simplified,E_X2_simplified,V_X_simplified,V_X_factored

# 2022 Q4(2)  2024.8.9

import numpy as np
import scipy.integrate as integrate

# パラメータの設定
gamma = 0.3  # 例として γ = 0.5 を設定
num_samples = 1000000  # サンプル数

# 確率密度関数 f(x) の定義
def pdf(x, gamma):
    return (1/gamma) * x**(-(1/gamma + 1))

# 理論的な期待値 E[X] の計算
def theoretical_expected_value(gamma):
    if gamma <= 0 or gamma >= 1:
        return np.inf
    else:
        result, _ = integrate.quad(lambda x: x * pdf(x, gamma), 1, np.inf)
        return result

# 理論的な分散 V[X] の計算
def theoretical_variance(gamma):
    if gamma <= 0 or gamma >= 0.5:
        return np.inf
    else:
        E_X = theoretical_expected_value(gamma)
        E_X2, _ = integrate.quad(lambda x: x**2 * pdf(x, gamma), 1, np.inf)
        return E_X2 - E_X**2

# シミュレーションによる確率変数 X のサンプリング
samples = (np.random.uniform(size=num_samples))**(-gamma)

# シミュレーションによる期待値 E[X] の推定
E_X_simulation = np.mean(samples)

# シミュレーションによる分散 V[X] の推定
V_X_simulation = np.var(samples)

# 理論値の計算
E_X_theoretical = theoretical_expected_value(gamma)
V_X_theoretical = theoretical_variance(gamma)

# 結果の表示
print(f"シミュレーションによる期待値 E[X]: {E_X_simulation}")
print(f"理論的な期待値 E[X]: {E_X_theoretical}")
print(f"シミュレーションによる分散 V[X]: {V_X_simulation}")
print(f"理論的な分散 V[X]: {V_X_theoretical}")

シミュレーションによる期待値 E[X]: 1.4284089260116615
理論的な期待値 E[X]: 1.4285714285714284
シミュレーションによる分散 V[X]: 0.4593987557711114
理論的な分散 V[X]: 0.45918367346938904

# 2022 Q4(1)  2024.8.8

import numpy as np
import matplotlib.pyplot as plt

def probability_density_function(x, gamma):
    if x <= 1:
        return 0
    else:
        return (1/gamma) * x**(-1/gamma-1)

# 配列のγ値に対するベクトル化された関数
def f_array(x, gammas):
    y_values = np.zeros((len(gammas), len(x)))
    for i, gamma in enumerate(gammas):
        y_values[i] = [probability_density_function(xi, gamma) for xi in x]
    return y_values

x_values = np.linspace(1.00001, 5, 500)
gammas = [0.5, 1, 2]
y_values_pdf_array = f_array(x_values, gammas)

plt.figure(figsize=(10, 7))

for i, gamma in enumerate(gammas):
    plt.plot(x_values, y_values_pdf_array[i], label=f"γ = {gamma}")

plt.title("確率密度関数 f(x) のグラフ")
plt.xlabel("x")
plt.ylabel("f(x)")
plt.grid(True)
plt.legend()
plt.show()

# 2022 Q4(1)  2024.8.8

import numpy as np
import matplotlib.pyplot as plt

# 累積分布関数を定義
def F(x, gamma):
    if x <= 1:
        return 0
    else:
        return 1 - x**(-1/gamma)

# 配列のγ値に対するベクトル化された関数
def F_array(x, gammas):
    y_values = np.zeros((len(gammas), len(x)))
    for i, gamma in enumerate(gammas):
        y_values[i] = [F(xi, gamma) for xi in x]
    return y_values

x_values = np.linspace(1, 5, 500)
gammas = [0.5, 1, 2]
y_values_array = F_array(x_values, gammas)

plt.figure(figsize=(10, 7))

for i, gamma in enumerate(gammas):
    plt.plot(x_values, y_values_array[i], label=f"γ = {gamma}")

plt.title("累積分布関数 F(x) のグラフ")
plt.xlabel("x")
plt.ylabel("F(x)")
plt.grid(True)
plt.legend()
plt.show()

# 2022 Q2(5)  2024.8.2

import numpy as np
from scipy.integrate import dblquad

# 確率密度関数
def pdf(u, v):
    return 1/3 if (u**2 - 2*u*v + v**2 <= 1) else 0

# 期待値 E[U] を計算
E_U, _ = dblquad(lambda v, u: u * pdf(u, v), -1, 1, lambda u: -1, lambda u: 1)
# 期待値 E[V] を計算
E_V, _ = dblquad(lambda v, u: v * pdf(u, v), -1, 1, lambda u: -1, lambda u: 1)

# 期待値 E[U^2] を計算
E_U2, _ = dblquad(lambda v, u: u**2 * pdf(u, v), -1, 1, lambda u: -1, lambda u: 1)
# 期待値 E[V^2] を計算
E_V2, _ = dblquad(lambda v, u: v**2 * pdf(u, v), -1, 1, lambda u: -1, lambda u: 1)

# 分散 Var(U) を計算
Var_U = E_U2 - E_U**2
# 分散 Var(V) を計算
Var_V = E_V2 - E_V**2

# 共分散 Cov(U, V) を計算
Cov_UV, _ = dblquad(lambda v, u: u*v * pdf(u, v), -1, 1, lambda u: -1, lambda u: 1)

# 相関係数 ρ(U, V) を計算
rho_UV = Cov_UV / np.sqrt(Var_U * Var_V)

E_U, E_V, Var_U, Var_V, rho_UV

(0.0, 0.0, 0.27777777931630593, 0.2777565638412289, 0.5000186020617616)

# 2022 Q2(5)  2024.8.2

import numpy as np

# モンテカルロ法の設定
num_points = 10000

# UとVの一様乱数を生成
U = np.random.uniform(-1, 1, num_points)
V = np.random.uniform(-1, 1, num_points)

# 条件 U^2 - 2UV + V^2 <= 1 を満たす点をフィルタリング
condition = U**2 - 2*U*V + V**2 <= 1
U_filtered = U[condition]
V_filtered = V[condition]

# フィルタリングされた点数
num_filtered_points = len(U_filtered)

# 期待値 E[U] と E[V] の計算
E_U = np.mean(U_filtered)
E_V = np.mean(V_filtered)

# 分散 Var(U) と Var(V) の計算
Var_U = np.var(U_filtered, ddof=1)
Var_V = np.var(V_filtered, ddof=1)

# 共分散 Cov(U, V) の計算
Cov_UV = np.cov(U_filtered, V_filtered, ddof=1)[0, 1]

# 相関係数 ρ(U, V) の計算
rho_UV = Cov_UV / np.sqrt(Var_U * Var_V)

E_U, E_V, Var_U, Var_V, rho_UV, num_filtered_points

(-0.0030089627151707833,
 -0.0028658643404494717,
 0.2780688116636205,
 0.2779098545694732,
 0.4835460989080292,
 7488)