import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# データ生成（正規分布 + 外れ値）
n = 200
data_normal = np.random.normal(loc=50, scale=10, size=n)
# 外れ値を追加
outliers = np.array([5, 10, 95, 100, 110])
data = np.concatenate([data_normal, outliers])

# IQRの計算
Q1 = np.percentile(data, 25)
Q2 = np.percentile(data, 50)
Q3 = np.percentile(data, 75)
IQR = Q3 - Q1

# 外れ値の検出
lower_fence = Q1 - 1.5 * IQR
upper_fence = Q3 + 1.5 * IQR
outlier_mask = (data < lower_fence) | (data > upper_fence)
n_outliers = np.sum(outlier_mask)

print(f"Q1 = {Q1:.2f}")
print(f"Q2 (中央値) = {Q2:.2f}")
print(f"Q3 = {Q3:.2f}")
print(f"IQR = {IQR:.2f}")
print(f"下限フェンス = {lower_fence:.2f}")
print(f"上限フェンス = {upper_fence:.2f}")
print(f"外れ値の数 = {n_outliers}")
print(f"外れ値: {data[outlier_mask]}")

# IQRと標準偏差の比較（ロバスト性のデモ）
# 汚染率を変えて比較
contamination_rates = np.linspace(0, 0.2, 20)
iqr_values = []
std_values = []

for rate in contamination_rates:
    n_contam = int(n * rate)
    clean_data = np.random.normal(50, 10, n)
    if n_contam > 0:
        contam_data = np.random.normal(150, 5, n_contam)  # 外れ値
        mixed_data = np.concatenate([clean_data, contam_data])
    else:
        mixed_data = clean_data
    iqr_values.append(np.percentile(mixed_data, 75) - np.percentile(mixed_data, 25))
    std_values.append(np.std(mixed_data))

# 可視化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 左上: ヒストグラムと四分位数
axes[0, 0].hist(data, bins=30, alpha=0.7, color='steelblue', edgecolor='black')
axes[0, 0].axvline(Q1, color='orange', linewidth=2, linestyle='--', label=f'Q1 = {Q1:.1f}')
axes[0, 0].axvline(Q2, color='red', linewidth=2, linestyle='-', label=f'Q2 = {Q2:.1f}')
axes[0, 0].axvline(Q3, color='green', linewidth=2, linestyle='--', label=f'Q3 = {Q3:.1f}')
axes[0, 0].axvline(lower_fence, color='purple', linewidth=1.5, linestyle=':',
                   label=f'Lower fence = {lower_fence:.1f}')
axes[0, 0].axvline(upper_fence, color='purple', linewidth=1.5, linestyle=':',
                   label=f'Upper fence = {upper_fence:.1f}')
axes[0, 0].set_xlabel('Value')
axes[0, 0].set_ylabel('Count')
axes[0, 0].set_title('Histogram with Quartiles and Fences')
axes[0, 0].legend(fontsize=8)

# 右上: 箱ひげ図
bp = axes[0, 1].boxplot(data, vert=True, patch_artist=True,
                        boxprops=dict(facecolor='lightblue', color='navy'),
                        medianprops=dict(color='red', linewidth=2),
                        flierprops=dict(marker='o', markerfacecolor='red', markersize=8))
axes[0, 1].set_ylabel('Value')
axes[0, 1].set_title('Box Plot')
# IQRの範囲を注釈
axes[0, 1].annotate('', xy=(1.15, Q3), xytext=(1.15, Q1),
                    arrowprops=dict(arrowstyle='<->', color='orange', lw=2))
axes[0, 1].text(1.2, (Q1+Q3)/2, f'IQR = {IQR:.1f}', fontsize=11, color='orange',
               va='center')

# 左下: 外れ値の検出結果
normal_mask = ~outlier_mask
axes[1, 0].scatter(np.arange(len(data))[normal_mask], data[normal_mask],
                  c='steelblue', s=15, alpha=0.5, label='Normal')
axes[1, 0].scatter(np.arange(len(data))[outlier_mask], data[outlier_mask],
                  c='red', s=50, marker='x', linewidths=2, label='Outlier')
axes[1, 0].axhline(upper_fence, color='purple', linestyle='--', alpha=0.5,
                   label=f'Upper fence ({upper_fence:.1f})')
axes[1, 0].axhline(lower_fence, color='purple', linestyle='--', alpha=0.5,
                   label=f'Lower fence ({lower_fence:.1f})')
axes[1, 0].set_xlabel('Index')
axes[1, 0].set_ylabel('Value')
axes[1, 0].set_title('Outlier Detection using IQR')
axes[1, 0].legend(fontsize=9)

# 右下: IQRとstdのロバスト性比較
axes[1, 1].plot(contamination_rates * 100, iqr_values / iqr_values[0], 'b-o',
               markersize=4, label='IQR (normalized)')
axes[1, 1].plot(contamination_rates * 100, std_values / std_values[0], 'r-o',
               markersize=4, label='Std (normalized)')
axes[1, 1].set_xlabel('Contamination Rate [%]')
axes[1, 1].set_ylabel('Normalized Value')
axes[1, 1].set_title('Robustness: IQR vs Standard Deviation')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()