import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
true_mean = 0
true_std = 1
sample_sizes = np.logspace(1, 3, 100).astype(int)  # 細かいサンプルサイズの範囲
n_simulations = 10000  # シミュレーション回数を増やす

# シミュレーション関数
def simulate_mean_estimator(size):
    samples = np.random.normal(true_mean, true_std, (n_simulations, size))
    return np.mean(samples, axis=1)

# 理論値とシミュレーション結果の計算
cr_lower_bounds = true_std**2 / sample_sizes
simulated_variances = [np.var(simulate_mean_estimator(size)) for size in sample_sizes]

# プロットの設定
plt.figure(figsize=(3, 2))
plt.plot(sample_sizes, cr_lower_bounds, 'r-', linewidth=0.3, label="Cramér-Rao Lower Bound")
plt.plot(sample_sizes, simulated_variances, 'b-', linewidth=0.3, alpha=0.7, label="Simulated Variance")
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Sample Size', fontsize=4)
plt.ylabel('Variance', fontsize=4)
plt.title('Cramér-Rao Lower Bound vs Simulated Variance\nfor Normal Distribution Mean Estimation', fontsize=4)
plt.legend(fontsize=4)
plt.grid(True, which="both", ls="-", alpha=0.2)

# 軸のtickを調整
plt.tick_params(axis='both', which='major', labelsize=4)

# 余白を調整
plt.tight_layout()

plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# パラメータ設定
true_mean = 5
true_std = 2
sample_sizes = [10, 50, 100, 500]
n_simulations = 10000

# 推定量の計算
def estimate_mean(data):
    return np.mean(data)

def estimate_median(data):
    return np.median(data)

# シミュレーション
results = {size: {'mean': [], 'median': []} for size in sample_sizes}

for size in sample_sizes:
    for _ in range(n_simulations):
        sample = np.random.normal(true_mean, true_std, size)
        results[size]['mean'].append(estimate_mean(sample))
        results[size]['median'].append(estimate_median(sample))

# プロットのスタイル設定
plt.style.use('seaborn')
plt.rcParams['font.size'] = 2
plt.rcParams['axes.labelsize'] = 2
plt.rcParams['axes.titlesize'] = 2
plt.rcParams['xtick.labelsize'] = 2
plt.rcParams['ytick.labelsize'] = 2
plt.rcParams['legend.fontsize'] = 2
plt.rcParams['figure.titlesize'] = 2

# プロット
fig, axs = plt.subplots(2, 2, figsize=(6, 6))
fig.suptitle('Sample Mean vs Sample Median as Estimators of Population Mean', fontsize=6)

for i, size in enumerate(sample_sizes):
    ax = axs[i//2, i%2]
    ax.hist(results[size]['mean'], bins=50, alpha=0.5, label='Sample Mean', color='skyblue', edgecolor='black')
    ax.hist(results[size]['median'], bins=50, alpha=0.5, label='Sample Median', color='lightgreen', edgecolor='black')
    ax.axvline(true_mean, color='red', linestyle='dashed', linewidth=0.5, label='True Mean')
    ax.set_title(f'Sample Size: {size}', fontsize=4)
    ax.set_xlabel('Estimated Value', fontsize=4)
    ax.set_ylabel('Frequency', fontsize=4)
    ax.legend(loc='upper left', fontsize=4)
    
    # x軸とy軸の範囲を設定
    ax.set_xlim(true_mean - 3*true_std/np.sqrt(size), true_mean + 3*true_std/np.sqrt(size))
    ax.set_ylim(0, n_simulations/8)  # ヒストグラムの最大高さを調整

    # グリッドを追加
    ax.grid(True, linestyle='--', alpha=0.7)

    # 軸のtickを調整
    ax.tick_params(axis='both', which='major', labelsize=4)
    
    
# サブプロット間の間隔を調整
plt.subplots_adjust(left=0.9, right=0.9, top=0.9, bottom=0.9, wspace=0.1, hspace=0.1)


plt.tight_layout(rect=[0, 0.1, 1, 0.85])  # タイトルのためのスペースを確保

# 統計量の計算と出力
print("\nEstimator Statistics:")
print("=====================")
for size in sample_sizes:
    print(f"\nSample Size: {size}")
    for estimator in ['mean', 'median']:
        estimates = results[size][estimator]
        bias = np.mean(estimates) - true_mean
        variance = np.var(estimates)
        mse = bias**2 + variance
        print(f"{estimator.capitalize()}:")
        print(f"  Bias: {bias:.4f}")
        print(f"  Variance: {variance:.4f}")
        print(f"  MSE: {mse:.4f}")

plt.show()