# 2018 Q3(5)  2024.10.17

from scipy.optimize import fsolve

# パラメータ設定
n = 8  # 試行回数
observed_mean = 3.5  # 仮定した標本平均

# モーメント法の方程式
def moment_eq(theta, observed_mean, n):
    return observed_mean - (n * theta) / (1 - (1 - theta)**n)

# 数値的に方程式を解く
theta_hat = fsolve(moment_eq, 0.1, args=(observed_mean, n))[0]

# 結果を表示
print(f"θ ハット: {theta_hat}")

θ ハット: 0.4328142334547646

# 2018 Q3(5)  2024.10.17

import numpy as np

# パラメータ設定
n = 8  # 試行回数
observed_mean = 3.5  # 仮定された標本平均
max_iterations = 100  # 最大反復回数
tolerance = 1e-6  # 収束判定の閾値

# 初期推定値
theta_hat = observed_mean / n

# 反復計算
for iteration in range(max_iterations):
    new_theta_hat = (observed_mean / n) * (1 - (1 - theta_hat) ** n)
    
    # 収束判定
    if abs(new_theta_hat - theta_hat) < tolerance:
        break
    
    theta_hat = new_theta_hat

# 結果の表示
print(f"θ ハット: {theta_hat}")

θ ハット: 0.432814320144486

# 2018 Q3(5)  2024.10.17

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
n = 8  # 試行回数
observed_mean = 3.5  # 仮定した標本平均
theta_values = np.linspace(0.01, 0.99, 100)  # θを0.01から0.99まで100点で変化させる

# f(θ) の定義
def f_theta(theta, n, observed_mean):
    return (n * theta) / (observed_mean * (1 - (1 - theta) ** n))

# f(θ) を計算
f_values = f_theta(theta_values, n, observed_mean)

# f(θ) = 1 となる交点を見つける
intersection_index = np.argmin(np.abs(f_values - 1))
theta_hat = theta_values[intersection_index]

# 交点の θ ハットの値を表示
print(f"θ ハット: {theta_hat}")

# グラフの描画
plt.plot(theta_values, f_values, label='f(θ)', color='blue')
plt.axhline(y=1, color='red', linestyle='--', label='y = 1')
plt.xlabel("θ", fontsize=12)
plt.ylabel("f(θ)", fontsize=12)
plt.title("f(θ) のグラフ (左辺 / 右辺の比)", fontsize=14)
plt.legend()
plt.grid(True)

# グラフの表示
plt.show()

θ ハット: 0.43565656565656563

# 2018 Q3(4)  2024.10.16

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
n = 8  # 試行回数
theta_values = np.linspace(0.01, 0.2, 50)  # θを0.01から0.2まで50点で変化

# 条件付き期待値 η(θ) と 条件なしの期待値 E[X] を計算
eta_values = n * theta_values / (1 - (1 - theta_values) ** n)  # η(θ)
EX_values = n * theta_values  # E[X]

# 二倍のE[X]をプロットするための値
double_EX_values = 2 * EX_values

# グラフの描画
plt.plot(theta_values, eta_values, label='条件付き期待値 η(θ)', color='blue', linestyle='-', marker='o', markersize=4)
plt.plot(theta_values, EX_values, label='条件なし期待値 E[X]', color='green', linestyle='--')
plt.plot(theta_values, double_EX_values, label='2倍の E[X]', color='red', linestyle=':')

# 二倍の関係になるθの位置を求め、プロット
intersection_index = np.argmin(np.abs(eta_values - double_EX_values))
intersection_theta = theta_values[intersection_index]
intersection_eta = eta_values[intersection_index]

plt.scatter(intersection_theta, intersection_eta, color='black', zorder=5, label=f'η(θ) = 2E[X] (θ ≈ {intersection_theta:.3f})')

# キャプションの表示
plt.xlabel("θ", fontsize=12)
plt.ylabel("期待値", fontsize=12)
plt.title("θ に対する η(θ) と E[X] の比較 (理論値)", fontsize=14)
plt.legend()
plt.grid(True)

# グラフの表示
plt.show()

# 2018 Q3(4)  2024.10.16

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
n = 8  # 試行回数
theta_values = np.linspace(0.01, 0.2, 50)  # θを0.01から0.2まで50点で変化
num_simulations = 100000  # シミュレーション回数

# シミュレーションによる条件付き期待値 η(θ) と 条件なし期待値 E[X] を計算
eta_values = []
EX_values = []

for theta in theta_values:
    # 二項分布のサンプリング
    binom_samples = np.random.binomial(n, theta, num_simulations)
    
    # 条件付きサンプル (X >= 1)
    samples_X_geq_1 = binom_samples[binom_samples >= 1]
    
    # 条件付き期待値 η(θ)
    eta = np.mean(samples_X_geq_1)
    eta_values.append(eta)
    
    # 条件なし期待値 E[X]
    EX = np.mean(binom_samples)
    EX_values.append(EX)

# 2倍のE[X]を計算
double_EX_values = 2 * np.array(EX_values)

# グラフの描画
plt.plot(theta_values, eta_values, label='条件付き期待値 η(θ)', color='blue', linestyle='-', marker='o', markersize=4)
plt.plot(theta_values, EX_values, label='条件なし期待値 E[X]', color='green', linestyle='--')
plt.plot(theta_values, double_EX_values, label='2倍の E[X]', color='red', linestyle=':')

# 二倍の関係になるθの位置を求め、プロット
intersection_index = np.argmin(np.abs(np.array(eta_values) - double_EX_values))
intersection_theta = theta_values[intersection_index]
intersection_eta = eta_values[intersection_index]

plt.scatter(intersection_theta, intersection_eta, color='black', zorder=5, label=f'η(θ) = 2E[X] (θ ≈ {intersection_theta:.3f})')

# キャプションの表示
plt.xlabel("θ", fontsize=12)
plt.ylabel("期待値", fontsize=12)
plt.title("θ に対する η(θ) と E[X] の比較 (シミュレーション)", fontsize=14)
plt.legend()
plt.grid(True)

# グラフの表示
plt.show()

# 2018 Q3(3)  2024.10.15

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
n = 10  # 試行回数
theta = 0.5  # 成功確率
num_simulations = 100000  # シミュレーション回数

# 二項分布のサンプリング
binom_samples = np.random.binomial(n, theta, num_simulations)

# 条件 X >= 1 を満たすサンプルを抽出
samples_X_geq_1 = binom_samples[binom_samples >= 1]

# 条件付き期待値と条件付き分散のシミュレーションによる計算
simulated_mean = np.mean(samples_X_geq_1)
simulated_variance = np.var(samples_X_geq_1)

# 条件なしの期待値と分散のシミュレーションによる計算
unconditional_simulated_mean = np.mean(binom_samples)
unconditional_simulated_variance = np.var(binom_samples)

# 理論値の計算
theoretical_mean = n * theta / (1 - (1 - theta) ** n)
theoretical_variance = (
    (n * theta * (1 - theta)) / (1 - (1 - theta) ** n)
    - (n ** 2 * theta ** 2 * (1 - theta) ** n) / (1 - (1 - theta) ** n) ** 2
)

# 条件なしの理論期待値と分散
unconditional_theoretical_mean = n * theta
unconditional_theoretical_variance = n * theta * (1 - theta)

# 結果を表示
print(f"シミュレーションによる条件付き期待値: {simulated_mean}")
print(f"理論上の条件付き期待値: {theoretical_mean}")
print(f"シミュレーションによる条件付き分散: {simulated_variance}")
print(f"理論上の条件付き分散: {theoretical_variance}")

print(f"シミュレーションによる条件なし期待値: {unconditional_simulated_mean}")
print(f"理論上の条件なし期待値: {unconditional_theoretical_mean}")
print(f"シミュレーションによる条件なし分散: {unconditional_simulated_variance}")
print(f"理論上の条件なし分散: {unconditional_theoretical_variance}")

# ヒストグラムの描画
plt.hist(samples_X_geq_1, bins=np.arange(1, n + 2) - 0.5, density=True, alpha=0.6, color='blue', edgecolor='black', label='条件付き (X >= 1)')
plt.hist(binom_samples, bins=np.arange(0, n + 2) - 0.5, density=True, alpha=0.4, color='green', edgecolor='black', label='条件なし')

# 理論値を折れ線グラフで表示
x_values = np.arange(0, n + 1)
binom_pmf = np.array([np.math.comb(n, x) * theta**x * (1 - theta)**(n - x) for x in x_values])
conditional_pmf = binom_pmf[1:] / (1 - binom_pmf[0])  # 条件付きのPMF

plt.plot(x_values[1:], conditional_pmf, marker='x', linestyle='-', color='red', label='条件付き理論値')
plt.plot(x_values, binom_pmf, marker='o', linestyle='-', color='purple', label='条件なし理論値')

# グラフの描画
plt.xlabel("X の値", fontsize=12)
plt.ylabel("確率密度", fontsize=12)
plt.title("条件付きと条件なしの二項分布のシミュレーション", fontsize=14)
plt.legend()
plt.grid(True)
plt.show()

シミュレーションによる条件付き期待値: 5.00641686604667
理論上の条件付き期待値: 5.004887585532747
シミュレーションによる条件付き分散: 2.4896127535349235
理論上の条件付き分散: 2.4779819766102995
シミュレーションによる条件なし期待値: 5.00106
理論上の条件なし期待値: 5.0
シミュレーションによる条件なし分散: 2.5137388763999997
理論上の条件なし分散: 2.5

シミュレーションによる条件付き期待値: 2.130760208977503
理論上の条件付き期待値: 2.1333333333333333
シミュレーションによる条件付き分散: 0.7789354600394489
理論上の条件付き分散: 0.7822222222222222
シミュレーションによる条件なし期待値: 1.99844
理論上の条件なし期待値: 2.0
シミュレーションによる条件なし分散: 0.9949975664
理論上の条件なし分散: 1.0

# 2018 Q3(3)  2024.10.15

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
theta = 0.5  # 成功確率（固定）
n_values = np.arange(1, 11, 1)  # nの値を1から10まで1刻みで変化させる
num_simulations = 100000  # シミュレーション回数

simulated_means = []  # 条件付き期待値のシミュレーション結果
theoretical_means = []  # 条件付き理論期待値
unconditional_means_sim = []  # 条件なしの期待値のシミュレーション結果
unconditional_means_theory = []  # 条件なしの理論期待値

# nを変化させてシミュレーション
for n in n_values:
    # 二項分布のサンプリング
    binom_samples = np.random.binomial(n, theta, num_simulations)
    samples_X_geq_1 = binom_samples[binom_samples >= 1]  # X >= 1 の条件に絞る
    
    # シミュレーションによる条件付き期待値
    simulated_mean = np.mean(samples_X_geq_1)
    simulated_means.append(simulated_mean)
    
    # 条件付き理論期待値
    theoretical_mean = n * theta / (1 - (1 - theta) ** n)
    theoretical_means.append(theoretical_mean)
    
    # 条件なしの期待値のシミュレーション結果
    unconditional_mean_sim = np.mean(binom_samples)
    unconditional_means_sim.append(unconditional_mean_sim)
    
    # 条件なしの理論期待値
    unconditional_mean_theory = n * theta
    unconditional_means_theory.append(unconditional_mean_theory)

# シミュレーション結果と理論値をプロット
plt.plot(n_values, simulated_means, label="条件付きシミュレーション", marker='o', linestyle='--', color='blue')
plt.plot(n_values, theoretical_means, label="条件付き理論値", marker='x', linestyle='-', color='red')
plt.plot(n_values, unconditional_means_sim, label="条件なしシミュレーション", marker='^', linestyle='--', color='purple')
plt.plot(n_values, unconditional_means_theory, label="条件なし理論値", marker='s', linestyle='-', color='green')

# グラフの描画
plt.xlabel("試行回数 n", fontsize=12)
plt.ylabel("期待値 E[X]", fontsize=12)
plt.title("n の変化に対する期待値の比較 (θ = 0.5)", fontsize=14)
plt.legend()
plt.grid(True)
plt.show()

# 2018 Q3(3)  2024.10.15

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
theta = 0.5  # 成功確率（固定）
n_values = np.arange(1, 11, 1)  # nの値を1から10まで1刻みで変化させる
num_simulations = 100000  # シミュレーション回数

simulated_variances = []  # 条件付き分散のシミュレーション結果
theoretical_variances = []  # 条件付き理論分散
unconditional_variances_sim = []  # 条件なしの分散のシミュレーション結果
unconditional_variances_theory = []  # 条件なしの理論分散

# nを変化させてシミュレーション
for n in n_values:
    # 二項分布のサンプリング
    binom_samples = np.random.binomial(n, theta, num_simulations)
    samples_X_geq_1 = binom_samples[binom_samples >= 1]  # X >= 1 の条件に絞る
    
    # シミュレーションによる条件付き分散
    simulated_variance = np.var(samples_X_geq_1)
    simulated_variances.append(simulated_variance)
    
    # 条件付き理論分散
    theoretical_variance = (
        (n * theta * (1 - theta)) / (1 - (1 - theta) ** n)
        - (n ** 2 * theta ** 2 * (1 - theta) ** n) / (1 - (1 - theta) ** n) ** 2
    )
    theoretical_variances.append(theoretical_variance)
    
    # 条件なしの分散のシミュレーション結果
    unconditional_variance_sim = np.var(binom_samples)
    unconditional_variances_sim.append(unconditional_variance_sim)
    
    # 条件なしの理論分散
    unconditional_variance_theory = n * theta * (1 - theta)
    unconditional_variances_theory.append(unconditional_variance_theory)

# シミュレーション結果と理論値をプロット
plt.plot(n_values, simulated_variances, label="条件付きシミュレーション", marker='o', linestyle='--', color='blue')
plt.plot(n_values, theoretical_variances, label="条件付き理論値", marker='x', linestyle='-', color='red')
plt.plot(n_values, unconditional_variances_sim, label="条件なしシミュレーション", marker='^', linestyle='--', color='purple')
plt.plot(n_values, unconditional_variances_theory, label="条件なし理論値", marker='s', linestyle='-', color='green')

# グラフの描画
plt.xlabel("試行回数 n", fontsize=12)
plt.ylabel("分散 V[X]", fontsize=12)
plt.title("n の変化に対する分散の比較 (θ = 0.5)", fontsize=14)
plt.legend()
plt.grid(True)
plt.show()

# 2018 Q3(2)  2024.10.12

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
n = 10  # 試行回数
theta = 0.5  # 成功確率
num_simulations = 10000  # シミュレーション回数

# 二項分布に従うランダム変数を生成
binom_samples = np.random.binomial(n, theta, num_simulations)

# X >= 1 の条件を満たすサンプルを抽出
samples_X_geq_1 = binom_samples[binom_samples >= 1]

# X = x の回数を数える
counts, _ = np.histogram(samples_X_geq_1, bins=np.arange(n+2))

# シミュレーションによる条件付き確率（X = 1 から n の範囲で計算）
simulated_probabilities = counts[1:] / np.sum(counts)

# 理論上の条件付き確率 h(x)
x_values = np.arange(1, n+1)
binom_pmf = np.array([np.math.comb(n, x) * theta**x * (1-theta)**(n-x) for x in x_values])
theoretical_probabilities = binom_pmf / (1 - (1 - theta)**n)

# シミュレーション結果をヒストグラムで表示
plt.bar(x_values, simulated_probabilities, color='blue', alpha=0.6, label="シミュレーション")
plt.plot(x_values, theoretical_probabilities, label="理論値", marker='x', linestyle='-', color='red')

# グラフの描画
plt.xlabel("X の値", fontsize=12)
plt.ylabel("条件付き確率 P(X = x | X >= 1)", fontsize=12)
plt.title(f"条件付き確率のシミュレーションと理論値の比較 (n={n}, θ={theta})", fontsize=14)
plt.legend()
plt.grid(True)
plt.show()