#2022 Q2(1)  2024.7.29

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

(u/2 + 1/2, v/2 + 1/2)

#2022 Q2(1)  2024.7.29

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

# 同時累積分布関数 F(u, v) を定義
def F(u, v):
    return (u * v + u + v + 1) / 4

# u, v の範囲を設定
u = np.linspace(-1, 1, 100)
v = np.linspace(-1, 1, 100)
u, v = np.meshgrid(u, v)
z = F(u, v)

# 3Dプロットの作成
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(u, v, z, cmap='viridis')

# グラフのラベル
ax.set_xlabel('U')
ax.set_ylabel('V')
ax.set_zlabel('F(U, V)')
ax.set_title('同時累積分布関数 F(U, V)')

plt.show()

#2022 Q1(4)  2024.7.28

# 定数
total_points = 16
PA = 3 / 4
PB = 3 / 4
PC = 3 / 4
PAB = 9 / 16
PAC = 9 / 16
PBC = 9 / 16

# 結果を格納するリスト
valid_probabilities = []

# 最小値と最大値を取るときの組み合わせを格納する変数
min_combination = None
max_combination = None
min_value = float('inf')
max_value = float('-inf')

# 各エリアの点数の範囲を狭める
for n1 in range(total_points // 3 + 1):
    for n2 in range((total_points - 3 * n1) // 3 + 1):
        for n3 in range(total_points - 3 * n1 - 3 * n2 + 1):
            n4 = total_points - 3 * n1 - 3 * n2 - n3
            if n4 < 0:
                continue

            calc_PA = (n1 + 2 * n2 + n4) / total_points
            calc_PB = (n1 + 2 * n2 + n4) / total_points
            calc_PC = (n1 + 2 * n2 + n4) / total_points
            calc_PAB = (n2 + n4) / total_points
            calc_PAC = (n2 + n4) / total_points
            calc_PBC = (n2 + n4) / total_points

            # 条件を満たすか確認
            if (abs(calc_PA - PA) < 1e-9 and 
                abs(calc_PB - PB) < 1e-9 and 
                abs(calc_PC - PC) < 1e-9 and 
                abs(calc_PAB - PAB) < 1e-9 and 
                abs(calc_PAC - PAC) < 1e-9 and 
                abs(calc_PBC - PBC) < 1e-9):
                
                current_PABC = n4 / total_points
                valid_probabilities.append(current_PABC)
                
                if current_PABC < min_value:
                    min_value = current_PABC
                    min_combination = (n1, n2, n3, n4)
                
                if current_PABC > max_value:
                    max_value = current_PABC
                    max_combination = (n1, n2, n3, n4)

# 最小値と最大値を表示
min_PABC = min_value if valid_probabilities else None
max_PABC = max_value if valid_probabilities else None

# 結果の表示
print(f"最小値: {min_PABC}")
print(f"最大値: {max_PABC}")
print(f"最小値を取るときのN1, N2, N3, N4: {min_combination}")
print(f"最大値を取るときのN1, N2, N3, N4: {max_combination}")

最小値: 0.375
最大値: 0.4375
最小値を取るときのN1, N2, N3, N4: (0, 3, 1, 6)
最大値を取るときのN1, N2, N3, N4: (1, 2, 0, 7)

### 2022 Q1(3)   2024.7.27

import itertools

def generate_combinations(n, k):
    base_list = [1] * k + [0] * (n - k)
    all_combinations = set(itertools.permutations(base_list))
    return list(all_combinations)

def calculate_intersection(A, B, C):
    return np.mean(np.array(A) & np.array(B) & np.array(C))

def find_min_max_intersection(n, k):
    all_combinations = generate_combinations(n, k)
    min_intersection = 1.0
    max_intersection = 0.0

    for A in all_combinations:
        for B in all_combinations:
            for C in all_combinations:
                intersection_ratio = calculate_intersection(A, B, C)
                if intersection_ratio < min_intersection:
                    min_intersection = intersection_ratio
                if intersection_ratio > max_intersection:
                    max_intersection = intersection_ratio

    return min_intersection, max_intersection

# パラメータを設定
n = 8
k = 6

# 最小値と最大値を計算
min_intersection, max_intersection = find_min_max_intersection(n, k)

min_intersection, max_intersection

(0.25, 0.75)

### 2022 Q1(1)   2024.7.26

import numpy as np

# シミュレーションの回数
num_simulations = 1000000

# 与えられた確率
prob = 3/4

# 事象 A, B, C をシミュレーション
A = np.random.rand(num_simulations) < prob
B = np.random.rand(num_simulations) < prob
C = np.random.rand(num_simulations) < prob

# 交差確率の計算
P_A_and_B = np.mean(A & B)
P_A_and_B_and_C = np.mean(A & B & C)

P_A_and_B, P_A_and_B_and_C

(0.56311, 0.422594)

#検算
9/16,27/64

(0.5625, 0.421875)

### 2022 Q1(2)   2024.7.26

import numpy as np

def simulate_overlap_with_fixed_true(correlation):
    P_A = 3/4
    P_B = 3/4
    num_elements = 1000
    num_true_A = int(P_A * num_elements)
    num_true_B = int(P_B * num_elements)

    # A と B の配列を生成し、3/4がtrueになるように設定
    A = np.zeros(num_elements, dtype=bool)
    A[:num_true_A] = True
    np.random.shuffle(A)

    B = np.zeros(num_elements, dtype=bool)
    B[:num_true_B] = True
    np.random.shuffle(B)

    if correlation == 'positive':
        # 正の相関を持たせるために、スワップを行う
        A1_B0_indices = np.where(A & ~B)[0]
        A0_B1_indices = np.where(~A & B)[0]

        min_len = min(len(A1_B0_indices), len(A0_B1_indices))
        for i in range(min_len):
            B[A1_B0_indices[i]], B[A0_B1_indices[i]] = B[A0_B1_indices[i]], B[A1_B0_indices[i]]

    elif correlation == 'negative':
        # 負の相関を持たせるために、スワップを行う
        A1_B1_indices = np.where(A & B)[0]
        A0_B0_indices = np.where(~A & ~B)[0]

        min_len = min(len(A1_B1_indices), len(A0_B0_indices))
        for i in range(min_len):
            B[A1_B1_indices[i]], B[A0_B0_indices[i]] = B[A0_B0_indices[i]], B[A1_B1_indices[i]]

    else:
        raise ValueError("Invalid correlation type")

    # 重なりの割合を計算
    overlap_ratio = np.mean(A & B)

    return overlap_ratio

# 正の相関の場合
overlap_positive = simulate_overlap_with_fixed_true('positive')

# 負の相関の場合
overlap_negative = simulate_overlap_with_fixed_true('negative')

overlap_positive, overlap_negative

(0.75, 0.5)

# 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