# 2021 Q5(2)  2024.9.6

import numpy as np

# 乱数で4つの要素を持つ基準ベクトルを生成 (整数で0以外の値のみ)
def generate_random_integer_vector_4(n=4):
    return np.random.randint(1, 10, size=n)  # 1から9の範囲で整数を生成

# 生成した基準ベクトルを使って、並行するベクトルを生成
def generate_parallel_matrix_4(v, scale_factors):
    L = np.array([scale * v for scale in scale_factors])
    return L

# ベクトル v に直交するベクトルを生成する関数
# グラム・シュミット法を適用して、基準ベクトルに直交するベクトルを作る
def generate_orthogonal_vectors_4(v):
    orthogonal_vectors = []
    for _ in range(5):
        candidate = np.random.randint(1, 10, size=v.shape)  # 新たな候補ベクトル
        # 直交条件：v と candidate の内積が 0 になるように調整
        adjustment = np.dot(candidate, v) / np.dot(v, v)
        orthogonal_vector = candidate - adjustment * v  # 直交ベクトルにする
        orthogonal_vectors.append(orthogonal_vector)
    return np.array(orthogonal_vectors)

# 基準ベクトルを生成
v_base = generate_random_integer_vector_4()

# 基準ベクトルと並行するベクトルを5つ作る
scale_factors_5 = [1, 2, 3, 4, 5]
L_5x4 = generate_parallel_matrix_4(v_base, scale_factors_5)

# 基準ベクトルに直交するベクトルを5つ作る
M_5x4 = generate_orthogonal_vectors_4(v_base)

# 行列 L, M の表示とその積の確認
L_5x4, M_5x4, np.dot(L_5x4, M_5x4.T)  # L, M, そしてその積を表示

(array([[ 7,  1,  6,  8],
        [14,  2, 12, 16],
        [21,  3, 18, 24],
        [28,  4, 24, 32],
        [35,  5, 30, 40]]),
 array([[-3.01333333,  8.42666667, -2.44      ,  3.41333333],
        [-1.68666667,  6.47333333,  3.84      , -2.21333333],
        [ 1.28666667,  4.32666667,  0.96      , -2.38666667],
        [ 1.14666667,  3.30666667,  2.84      , -3.54666667],
        [-2.73333333,  4.46666667, -1.2       ,  2.73333333]]),
 array([[-7.10542736e-15,  7.10542736e-15,  0.00000000e+00,
         -3.55271368e-15,  0.00000000e+00],
        [-1.42108547e-14,  1.42108547e-14,  0.00000000e+00,
         -7.10542736e-15,  0.00000000e+00],
        [-2.13162821e-14,  2.13162821e-14,  7.10542736e-15,
          0.00000000e+00, -7.10542736e-15],
        [-2.84217094e-14,  2.84217094e-14,  0.00000000e+00,
         -1.42108547e-14,  0.00000000e+00],
        [-3.55271368e-14,  2.84217094e-14,  0.00000000e+00,
         -2.84217094e-14, -7.10542736e-15]]))

# 2021 Q5(2)  2024.9.6

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# l, m, n を定義
l = 5  # Lの行数
m = 5  # Mの行数
n = 4  # LとMの列数 (xの次元)

# L と M を定義 (与えられた行列を使用)
L = np.array([[7, 1, 6, 8],
              [14, 2, 12, 16],
              [21, 3, 18, 24],
              [28, 4, 24, 32],
              [35, 5, 30, 40]])

M = np.array([[-3.01333333, 8.42666667, -2.44, 3.41333333],
              [-1.68666667, 6.47333333, 3.84, -2.21333333],
              [1.28666667, 4.32666667, 0.96, -2.38666667],
              [1.14666667, 3.30666667, 2.84, -3.54666667],
              [-2.73333333, 4.46666667, -1.2, 2.73333333]])

# N(0, 1) に従う乱数ベクトル x を生成 (サンプル数を10000000にする)
sample_size = 10000000
x = np.random.randn(n, sample_size)

# Lx と Mx を計算
Lx = np.dot(L, x)  # lxサンプル数の行列
Mx = np.dot(M, x)  # mxサンプル数の行列

# 共分散行列 Cov(Lx, Mx) を計算
cov_matrix = np.zeros((l, m))
for i in range(l):
    for j in range(m):
        cov_matrix[i, j] = np.cov(Lx[i, :], Mx[j, :])[0, 1]

# L * M^T の積を再計算
LM_T = np.dot(L, M.T)

# ヒートマップを描画
plt.figure(figsize=(12, 6))

# 共分散行列のヒートマップ
plt.subplot(1, 2, 1)
sns.heatmap(cov_matrix, cmap="coolwarm", annot=True, fmt=".2f", center=0, vmin=-1, vmax=1)
plt.title(f'共分散行列 Cov(Lx, Mx)')

# L * M^T のヒートマップ
plt.subplot(1, 2, 2)
sns.heatmap(LM_T, cmap="coolwarm", annot=True, fmt=".2f", center=0, vmin=-1, vmax=1)
plt.title(f'行列 L * M^T')

plt.tight_layout()
plt.show()

# 2021 Q5(1)  2024.9.4

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# 標準正規分布からサンプルを生成
n = 1000  # サンプル数
d = 5  # 次元数
X = np.random.randn(d, n)

# 4x5の線形変換行列 L を定義
L = np.array([[9, 7, 8, 6, 5],
              [6, 5, 4, 3, 2],
              [3, 2, 1, 9, 8],
              [8, 7, 6, 5, 0]])

# 線形変換を適用して y を計算
Y = np.dot(L, X)

# y のサンプルから分散共分散行列を推定
V_Y_estimated = np.cov(Y)

# 理論的な L・L^T と比較
V_Y_theoretical = np.dot(L, L.T)

# ヒートマップのプロット
plt.figure(figsize=(12, 6))

# 推定された分散共分散行列のヒートマップ
plt.subplot(1, 2, 1)
sns.heatmap(V_Y_estimated, annot=True, cmap="coolwarm", cbar=True)
plt.title("推定された分散共分散行列")

# 理論的な分散共分散行列のヒートマップ
plt.subplot(1, 2, 2)
sns.heatmap(V_Y_theoretical, annot=True, cmap="coolwarm", cbar=True)
plt.title("理論的な分散共分散行列")

plt.show()

# 2021 Q4(5)  2024.9.3

import numpy as np

def find_ai_sequence(n, tolerance=1e-6):
    while True:
        # ランダムに n 個の値を生成
        a = np.random.uniform(-1, 1, n)
        
        # 正の数と負の数の個数を調整
        neg_count = (n + 1) // 2  # 負の数を1つ多くする
        pos_count = n - neg_count
        
        a[:pos_count] = np.abs(a[:pos_count])
        a[pos_count:] = -np.abs(a[pos_count:])
        
        # 合計が0になるように調整
        a -= np.mean(a)
        
        # 条件をチェック
        sum_ai = np.sum(a)
        sum_ai_cubed = np.sum(a**3)
        
        if np.abs(sum_ai) < tolerance and np.abs(sum_ai_cubed - 1) < tolerance:
            return a

# パラメータの設定
n = 10  # サンプル数

# 条件を満たす a_i の数列を見つける
a = find_ai_sequence(n)

# 計算を確認
print(f"a_i の合計 (≒0): {a.sum()}")
print(f"a_i^3 の合計 (≒1): {np.sum(a ** 3)}")

# 条件を満たした a_i を表示
print("a_i の数列:")
print(a)

a_i の合計 (≒0): 6.661338147750939e-16
a_i^3 の合計 (≒1): 1.0000001485058854
a_i の数列:
[ 0.13279542  1.05885499  0.24951311  0.83492921  0.13800083 -0.71178855
 -0.38213634 -0.18394766 -0.55820111 -0.5780199 ]

import numpy as np
import matplotlib.pyplot as plt

# すでに見つかった条件を満たす a_i を使用
a = np.array([ 0.13279542,  1.05885499,  0.24951311,  0.83492921,  0.13800083, -0.71178855,
              -0.38213634, -0.18394766, -0.55820111, -0.5780199 ])

# パラメータの設定
n = len(a)  # サンプル数（a_i の数と一致させる）
scale = 1.0  # 指数分布のスケールパラメータ
iterations = 10000  # シミュレーション回数

# シミュレーション結果を保存するリスト
tau_hat_values = []

# シミュレーションを実行
for _ in range(iterations):
    # 指数分布に従う乱数を生成
    X = np.random.exponential(scale=scale, size=n)
    
    # \hat{\tau} の計算
    tau_hat = np.sum(a * X) ** 3
    tau_hat_values.append(tau_hat)

# 指数分布の三次中心モーメントを計算
tau = 2 * scale**3  # 指数分布における三次中心モーメント

# シミュレーション結果の平均
simulated_mean_tau_hat = np.mean(tau_hat_values)

# 結果を表示
print(f"理論値の τ: {tau}")
print(f"シミュレーションによる τ の推定値: {simulated_mean_tau_hat}")

# 視覚化
plt.figure(figsize=(10, 6))
plt.hist(tau_hat_values, bins=50, alpha=0.75, label='シミュレーション結果')
plt.axvline(simulated_mean_tau_hat, color='red', linestyle='dashed', linewidth=2, label=f'シミュレーション平均 = {simulated_mean_tau_hat:.2f}')
plt.axvline(tau, color='blue', linestyle='dashed', linewidth=2, label=f'理論値 = {tau:.2f}')
plt.xlabel('値')
plt.ylabel('頻度')
plt.title(r"$\hat{\tau} = \left(\sum_{i=1}^{n} a_i X_i\right)^3$ の分布 (指数分布)")
plt.legend()
plt.grid(True)
plt.show()

理論値の τ: 2.0
シミュレーションによる τ の推定値: 2.0555706049874143

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
n = 10  # サンプル数 (必要に応じて変更可能)
scale = 1.0  # 指数分布のスケールパラメータ（1/λ）
iterations = 10000  # シミュレーション回数

# シミュレーション結果を保存するリスト
sum_X_i_minus_X_bar_cubed_values = []

# シミュレーションを実行
for _ in range(iterations):
    # 指数分布に従う乱数を生成
    X = np.random.exponential(scale=scale, size=n)
    
    # 標本平均を計算
    X_bar = np.mean(X)
    
    # (X_i - X_bar)^3 の合計を計算
    sum_X_i_minus_X_bar_cubed = np.sum((X - X_bar)**3)
    sum_X_i_minus_X_bar_cubed_values.append(sum_X_i_minus_X_bar_cubed)

# 指数分布の三次中心モーメントを計算
tau = 2 * scale**3  # 指数分布における三次中心モーメント

# 理論値の計算
theoretical_value = (n - 1) * (n - 2) * tau / n

# シミュレーション結果の平均
simulated_mean = np.mean(sum_X_i_minus_X_bar_cubed_values)

# 結果を表示
print(f"理論値: {theoretical_value}")
print(f"シミュレーションによる平均: {simulated_mean}")

# 視覚化
plt.figure(figsize=(10, 6))
plt.hist(sum_X_i_minus_X_bar_cubed_values, bins=50, alpha=0.75, label='シミュレーション結果')
plt.axvline(simulated_mean, color='red', linestyle='dashed', linewidth=2, label=f'シミュレーション平均 = {simulated_mean:.2f}')
plt.axvline(theoretical_value, color='blue', linestyle='dashed', linewidth=2, label=f'理論値 = {theoretical_value:.2f}')
plt.xlabel('値')
plt.ylabel('頻度')
plt.title(r"$\sum_{i=1}^{n} (X_i - \overline{X})^3$ の分布 (指数分布)")
plt.legend()
plt.grid(True)
plt.show()

理論値: 14.4
シミュレーションによる平均: 14.499748799877718

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
n = 1000  # サンプル数
mu, sigma = 0, 1  # 正規分布の平均と標準偏差
iterations = 10000  # シミュレーション回数

# シミュレーション結果を保存するリスト
Y_i_cubed_values = []
Y_i_squared_Y_j_values = []
Y_i_Y_j_Y_k_values = []

# シミュレーションを実行
for _ in range(iterations):
    X_i = np.random.normal(mu, sigma, size=n)
    X_j = np.random.normal(mu, sigma, size=n)
    X_k = np.random.normal(mu, sigma, size=n)
    
    Y_i = X_i - mu
    Y_j = X_j - mu
    Y_k = X_k - mu
    
    Y_i_cubed_values.append(np.mean(Y_i**3))
    Y_i_squared_Y_j_values.append(np.mean(Y_i**2 * Y_j))
    Y_i_Y_j_Y_k_values.append(np.mean(Y_i * Y_j * Y_k))

# グラフを横に並べて表示
fig, axes = plt.subplots(1, 3, figsize=(18, 6))

# (1) E[Y_i^3] の分布
axes[0].hist(Y_i_cubed_values, bins=50, alpha=0.75)
axes[0].axvline(np.mean(Y_i_cubed_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[0].set_title(r"正規分布における $E[Y_i^3]$ の分布")
axes[0].set_xlabel('値')
axes[0].set_ylabel('頻度')
axes[0].legend()

# (2) E[Y_i^2 Y_j] の分布
axes[1].hist(Y_i_squared_Y_j_values, bins=50, alpha=0.75)
axes[1].axvline(np.mean(Y_i_squared_Y_j_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[1].set_title(r"正規分布における $E[Y_i^2 Y_j]$ の分布")
axes[1].set_xlabel('値')
axes[1].set_ylabel('頻度')
axes[1].legend()

# (3) E[Y_i Y_j Y_k] の分布
axes[2].hist(Y_i_Y_j_Y_k_values, bins=50, alpha=0.75)
axes[2].axvline(np.mean(Y_i_Y_j_Y_k_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[2].set_title(r"正規分布における $E[Y_i Y_j Y_k]$ の分布")
axes[2].set_xlabel('値')
axes[2].set_ylabel('頻度')
axes[2].legend()

# グラフを表示
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
n = 1000  # サンプル数
scale = 1.0  # 指数分布のスケールパラメータ（1/λ）
iterations = 10000  # シミュレーション回数

# シミュレーション結果を保存するリスト
Y_i_cubed_values = []
Y_i_squared_Y_j_values = []
Y_i_Y_j_Y_k_values = []

# シミュレーションを実行
for _ in range(iterations):
    X_i = np.random.exponential(scale=scale, size=n)
    X_j = np.random.exponential(scale=scale, size=n)
    X_k = np.random.exponential(scale=scale, size=n)
    
    Y_i = X_i - np.mean(X_i)
    Y_j = X_j - np.mean(X_j)
    Y_k = X_k - np.mean(X_k)
    
    Y_i_cubed_values.append(np.mean(Y_i**3))
    Y_i_squared_Y_j_values.append(np.mean(Y_i**2 * Y_j))
    Y_i_Y_j_Y_k_values.append(np.mean(Y_i * Y_j * Y_k))

# グラフを横に並べて表示
fig, axes = plt.subplots(1, 3, figsize=(18, 6))

# (1) E[Y_i^3] の分布
axes[0].hist(Y_i_cubed_values, bins=50, alpha=0.75)
axes[0].axvline(np.mean(Y_i_cubed_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[0].set_title(r"指数分布における $E[Y_i^3]$ の分布")
axes[0].set_xlabel('値')
axes[0].set_ylabel('頻度')
axes[0].legend()

# (2) E[Y_i^2 Y_j] の分布
axes[1].hist(Y_i_squared_Y_j_values, bins=50, alpha=0.75)
axes[1].axvline(np.mean(Y_i_squared_Y_j_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[1].set_title(r"指数分布における $E[Y_i^2 Y_j]$ の分布")
axes[1].set_xlabel('値')
axes[1].set_ylabel('頻度')
axes[1].legend()

# (3) E[Y_i Y_j Y_k] の分布
axes[2].hist(Y_i_Y_j_Y_k_values, bins=50, alpha=0.75)
axes[2].axvline(np.mean(Y_i_Y_j_Y_k_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[2].set_title(r"指数分布における $E[Y_i Y_j Y_k]$ の分布")
axes[2].set_xlabel('値')
axes[2].set_ylabel('頻度')
axes[2].legend()

# グラフを表示
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
n = 1000  # サンプル数
shape, scale = 2.0, 2.0  # ガンマ分布の形状パラメータとスケールパラメータ
iterations = 10000  # シミュレーション回数

# シミュレーション結果を保存するリスト
Y_i_cubed_values = []
Y_i_squared_Y_j_values = []
Y_i_Y_j_Y_k_values = []

# シミュレーションを実行
for _ in range(iterations):
    X_i = np.random.gamma(shape, scale, size=n)
    X_j = np.random.gamma(shape, scale, size=n)
    X_k = np.random.gamma(shape, scale, size=n)
    
    Y_i = X_i - np.mean(X_i)
    Y_j = X_j - np.mean(X_j)
    Y_k = X_k - np.mean(X_k)
    
    Y_i_cubed_values.append(np.mean(Y_i**3))
    Y_i_squared_Y_j_values.append(np.mean(Y_i**2 * Y_j))
    Y_i_Y_j_Y_k_values.append(np.mean(Y_i * Y_j * Y_k))

# グラフを横に並べて表示
fig, axes = plt.subplots(1, 3, figsize=(18, 6))

# (1) E[Y_i^3] の分布
axes[0].hist(Y_i_cubed_values, bins=50, alpha=0.75)
axes[0].axvline(np.mean(Y_i_cubed_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[0].set_title(r"ガンマ分布における $E[Y_i^3]$ の分布")
axes[0].set_xlabel('値')
axes[0].set_ylabel('頻度')
axes[0].legend()

# (2) E[Y_i^2 Y_j] の分布
axes[1].hist(Y_i_squared_Y_j_values, bins=50, alpha=0.75)
axes[1].axvline(np.mean(Y_i_squared_Y_j_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[1].set_title(r"ガンマ分布における $E[Y_i^2 Y_j]$ の分布")
axes[1].set_xlabel('値')
axes[1].set_ylabel('頻度')
axes[1].legend()

# (3) E[Y_i Y_j Y_k] の分布
axes[2].hist(Y_i_Y_j_Y_k_values, bins=50, alpha=0.75)
axes[2].axvline(np.mean(Y_i_Y_j_Y_k_values), color='red', linestyle='dashed', linewidth=2, label='平均')
axes[2].set_title(r"ガンマ分布における $E[Y_i Y_j Y_k]$ の分布")
axes[2].set_xlabel('値')
axes[2].set_ylabel('頻度')
axes[2].legend()

# グラフを表示
plt.tight_layout()
plt.show()

# 2021 Q4(1)  2024.8.31

import numpy as np

# パラメータの設定
n = 100  # サンプル数
mu = 0  # 平均
sigma = 1  # 標準偏差
iterations = 10000  # シミュレーション回数

# 各シミュレーションで得た標本平均を格納するリスト
sample_means = []

for _ in range(iterations):
    # 正規分布に従う乱数を生成
    X = np.random.normal(loc=mu, scale=sigma, size=n)
    
    # 標本平均を計算
    X_bar = np.mean(X)
    sample_means.append(X_bar)

# 標本平均の分散を計算
estimated_variance_of_sample_mean = np.var(sample_means)

# 理論的な分散
theoretical_variance = sigma**2 / n

print(f"理論的な分散: {theoretical_variance}")
print(f"推定された標本平均の分散: {estimated_variance_of_sample_mean}")

理論的な分散: 0.01
推定された標本平均の分散: 0.009838287336429174

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
n = 100  # サンプル数
mu = 0  # 平均
sigma = 1  # 標準偏差
iterations = 10000  # シミュレーション回数

# 分散の推定値を格納するリスト
estimated_variances = []

for _ in range(iterations):
    # 正規分布に従う乱数を生成
    X = np.random.normal(loc=mu, scale=sigma, size=n)
    
    # 標本平均の分散を計算
    estimated_variance = np.var(X) / n
    estimated_variances.append(estimated_variance)

# 理論的な分散
theoretical_variance = sigma**2 / n

# 視覚化
plt.figure(figsize=(10, 6))
plt.hist(estimated_variances, bins=100, alpha=0.75, label='推定された分散', range=(0.005, 0.015))
plt.axvline(theoretical_variance, color='red', linestyle='dashed', linewidth=2, label=f'理論的な分散 = {theoretical_variance}')
plt.xlabel('分散')
plt.ylabel('頻度')
plt.title(f'{iterations}回のシミュレーションにおける推定された分散の分布\n(サンプルサイズ = {n})')
plt.legend()
plt.grid(True)
plt.show()

# 2021 Q4(1)  2024.8.31

import numpy as np

# パラメータの設定
n_simulations = 10000  # シミュレーション回数
mu = 0  # 平均
sigma = 1  # 標準偏差

# X1, X2を正規分布から生成
X1 = np.random.normal(loc=mu, scale=sigma, size=n_simulations)
X2 = np.random.normal(loc=mu, scale=sigma, size=n_simulations)

# (X1 - X2)^3を計算
expression = (X1 - X2)**3

# 期待値を計算
expected_value = np.mean(expression)

print(f"(X1 - X2)^3 の期待値: {expected_value}")

(X1 - X2)^3 の期待値: 0.044962030052379295

import numpy as np
import matplotlib.pyplot as plt

# パラメータの設定
n = 10000  # シミュレーション回数
mu = 0  # 平均
sigma = 1  # 標準偏差

# (X1 - X2)^3 の推定値を格納するリスト
expected_values = []

for _ in range(n):
    # 正規分布に従う X1, X2 を生成
    X1 = np.random.normal(loc=mu, scale=sigma)
    X2 = np.random.normal(loc=mu, scale=sigma)
    
    # (X1 - X2)^3 を計算
    value = (X1 - X2)**3
    expected_values.append(value)

# 理論的な期待値
theoretical_expectation = 0

# 視覚化
plt.figure(figsize=(10, 6))
plt.hist(expected_values, bins=100, alpha=0.75, label='(X1 - X2)^3 の推定値', range=(-30, 30))
plt.axvline(theoretical_expectation, color='red', linestyle='dashed', linewidth=2, label=f'理論的な期待値 = {theoretical_expectation}')
plt.xlabel('値')
plt.ylabel('頻度')
plt.title(f'{n}回のシミュレーションにおける (X1 - X2)^3 の期待値の分布')
plt.legend()
plt.grid(True)
plt.show()