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

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