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

#2022 Q2(4) 2024.8.1

import numpy as np
import scipy.integrate as integrate

# 積分の被積分関数を定義
def integrand(v, u):
    return 1/4 if (u**2 - 2*u*v + v**2 <= 1) else 0

# 積分の範囲を定義する関数
def bounds_u():
    return -1, 1

def bounds_v(u):
    return -1, 1

# 確率 P(U^2 - 2UV + V^2 <= 1) を計算
probability, _ = integrate.nquad(integrand, [bounds_v, bounds_u])

probability

0.7500047152195971

#2022 Q2(4) 2024.8.1

import matplotlib.pyplot as plt
import numpy as np

# UとVの範囲を設定
u = np.linspace(-1, 1, 400)
v = np.linspace(-1, 1, 400)
U, V = np.meshgrid(u, v)

# 条件 U^2 - 2UV + V^2 <= 1 を満たす領域を計算
condition = (U**2 - 2*U*V + V**2 <= 1)

# グラフの描画
plt.figure(figsize=(5, 5))
plt.contourf(U, V, condition, levels=[0.5, 1], colors=['lightblue'], alpha=0.5)
plt.title(r'条件 $U^2 - 2UV + V^2 \leq 1$')
plt.xlabel('U')
plt.ylabel('V')
plt.axhline(0, color='gray', lw=0.5)
plt.axvline(0, color='gray', lw=0.5)
plt.grid(True)

# 条件を満たす境界を描画
plt.contour(U, V, condition, levels=[0.5], colors=['blue'])

plt.show()

#2022 Q2(2) 2024.7.30

from sympy import symbols, integrate

# 変数の定義
u, v = symbols('u v')

# 同時累積分布関数の定義
F_uv = (u*v + u + v + 1) / 4

# 周辺累積分布関数 F1(u) と F2(v) の計算
F1_u = integrate(F_uv, (v, -1, 1))
F2_v = integrate(F_uv, (u, -1, 1))

# 周辺確率密度関数 f1(u) と f2(v) の計算
f1_u = F1_u.diff(u)
f2_v = F2_v.diff(v)

# 同時確率密度関数 f(u,v) を求める
f_uv = F_uv.diff(u).diff(v)

# f(u,v) を f1(u) と f2(v) の積として表す
product_f1_f2 = f1_u * f2_v

# 結果の出力
results = {
    'F1_u': F1_u,
    'F2_v': F2_v,
    'f1_u': f1_u,
    'f2_v': f2_v,
    'f_uv': f_uv,
    'product_f1_f2': product_f1_f2
}

results

{'F1_u': u/2 + 1/2,
 'F2_v': v/2 + 1/2,
 'f1_u': 1/2,
 'f2_v': 1/2,
 'f_uv': 1/4,
 'product_f1_f2': 1/4}

#2022 Q2(3) 2024.7.31

from sympy import symbols, pi
import numpy as np
from scipy.integrate import quad

# 符号変数の定義
u, v = symbols('u v')

# f(u,v) の定義
f_uv_value = 1/4

# 条件を満たす範囲でのf(u,v)（数値積分用）
def integrand(v, u):
    if u**2 + v**2 <= 1:
        return f_uv_value
    else:
        return 0

# uに関しての積分
def integrate_u(u):
    return quad(integrand, -1, 1, args=(u))[0]

# 数値積分による評価
probability_numeric = quad(integrate_u, -1, 1)[0]

# 手計算による値 (π/4) の計算
manual_calculation_value = pi / 4

# 結果の出力
print(f"数値積分による確率: {probability_numeric}\n")
print(f"手計算による値 (π/4): {manual_calculation_value.evalf()}\n")

数値積分による確率: 0.7855319594076084

手計算による値 (π/4): 0.785398163397448

#2022 Q2(3) 2024.7.31

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# u, vの範囲と値を定義
u = np.linspace(-1, 1, 100)
v = np.linspace(-1, 1, 100)
U, V = np.meshgrid(u, v)

# f(u,v)の値を計算
f_uv_value = 1/4
F = np.where(U**2 + V**2 <= 1, f_uv_value, 0)

# 3Dプロットを作成
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(U, V, F, cmap='viridis')
ax.set_xlabel('U')
ax.set_ylabel('V')
ax.set_zlabel('f(u,v)')
ax.set_zlim(0, 1)  # z軸の範囲を設定
ax.set_title('領域 U^2 + V^2 <= 1 における f(u,v) の3Dプロット')
plt.show()