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

def ks_test_one_sample(data, cdf_func):
    """
    一標本KS検定のスクラッチ実装

    Parameters
    ----------
    data : array-like  観測データ
    cdf_func : callable  帰無仮説のCDF

    Returns
    -------
    D : float  KS統計量
    D_plus, D_minus : float  片側統計量
    max_idx : int  最大乖離の位置
    """
    x = np.sort(data)
    n = len(x)

    # F_0の値
    F0 = cdf_func(x)

    # D+とD-
    D_plus_vals = np.arange(1, n + 1) / n - F0
    D_minus_vals = F0 - np.arange(0, n) / n

    D_plus = np.max(D_plus_vals)
    D_minus = np.max(D_minus_vals)
    D = max(D_plus, D_minus)

    if D_plus >= D_minus:
        max_idx = np.argmax(D_plus_vals)
    else:
        max_idx = np.argmax(D_minus_vals)

    return D, D_plus, D_minus, max_idx

# --- テストデータ ---
np.random.seed(42)
n = 30
data = np.random.normal(loc=0.3, scale=1.0, size=n)  # 少しずれた正規分布

# H0: N(0, 1) に従う
D, Dp, Dm, idx = ks_test_one_sample(data, stats.norm.cdf)

# SciPyとの比較
D_scipy, p_scipy = stats.kstest(data, "norm")

print(f"KS統計量:  D = {D:.4f} (SciPy: {D_scipy:.4f})")
print(f"D+ = {Dp:.4f}, D- = {Dm:.4f}")
print(f"p値 (SciPy): {p_scipy:.6f}")

# --- 可視化 ---
fig, ax = plt.subplots(figsize=(10, 6))

# ソートされたデータ
x_sorted = np.sort(data)
n = len(x_sorted)

# 経験分布関数（階段関数）
for i in range(n):
    x_left = x_sorted[i - 1] if i > 0 else x_sorted[0] - 1
    x_right = x_sorted[i]
    ax.plot([x_left, x_right], [i / n, i / n], "b-", linewidth=1.5)
    ax.plot([x_right, x_right], [i / n, (i + 1) / n], "b-", linewidth=1.5)
# 最後のステップ
ax.plot([x_sorted[-1], x_sorted[-1] + 1.5], [1, 1], "b-", linewidth=1.5,
        label=r"$F_n(x)$ (EDF)")

# 理論分布
x_theo = np.linspace(x_sorted[0] - 1, x_sorted[-1] + 1.5, 500)
ax.plot(x_theo, stats.norm.cdf(x_theo), "r-", linewidth=2,
        label=r"$F_0(x) = \Phi(x)$ (Standard Normal)")

# 最大乖離の表示
x_max = x_sorted[idx]
F_n_val = (idx + 1) / n if Dp >= Dm else idx / n
F_0_val = stats.norm.cdf(x_max)
ax.annotate("", xy=(x_max, F_n_val), xytext=(x_max, F_0_val),
            arrowprops=dict(arrowstyle="<->", color="green", lw=2))
ax.text(x_max + 0.15, (F_n_val + F_0_val) / 2, f"D = {D:.3f}",
        fontsize=12, color="green", fontweight="bold")

ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("Cumulative probability", fontsize=12)
ax.set_title(f"KS test: EDF vs Standard Normal (n={n})", fontsize=13)
ax.legend(fontsize=11, loc="upper left")
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("ks_test_one_sample.png", dpi=150, bbox_inches="tight")
plt.show()

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

np.random.seed(42)

def kolmogorov_cdf(t, n_terms=100):
    """コルモゴロフ分布のCDF"""
    result = 0.0
    for j in range(1, n_terms + 1):
        result += (-1)**(j - 1) * np.exp(-2 * j**2 * t**2)
    return 1 - 2 * result

n_list = [10, 30, 100, 500]
n_sim = 10000

fig, axes = plt.subplots(2, 2, figsize=(12, 10))

for ax, n in zip(axes.flat, n_list):
    # シミュレーション: 帰無仮説のもとでKS統計量を計算
    D_values = np.zeros(n_sim)
    for i in range(n_sim):
        data = np.random.uniform(0, 1, n)  # U(0,1)から生成
        data_sorted = np.sort(data)
        F0 = data_sorted  # U(0,1)のCDF = x
        D_plus = np.max(np.arange(1, n + 1) / n - F0)
        D_minus = np.max(F0 - np.arange(0, n) / n)
        D_values[i] = max(D_plus, D_minus)

    # sqrt(n) * D の分布
    scaled_D = np.sqrt(n) * D_values

    ax.hist(scaled_D, bins=50, density=True, alpha=0.6, color="steelblue",
            edgecolor="white", label=f"Simulation (n={n})")

    # コルモゴロフ分布の密度関数（数値微分で近似）
    t = np.linspace(0.2, 2.5, 200)
    K_cdf = np.array([kolmogorov_cdf(ti) for ti in t])
    K_pdf = np.gradient(K_cdf, t)
    ax.plot(t, K_pdf, "r-", linewidth=2, label="Kolmogorov dist.")

    ax.set_xlabel(r"$\sqrt{n} \cdot D_n$", fontsize=11)
    ax.set_ylabel("Density", fontsize=11)
    ax.set_title(f"n = {n}", fontsize=12)
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)
    ax.set_xlim(0, 2.5)

plt.suptitle(r"Convergence of $\sqrt{n}D_n$ to Kolmogorov distribution",
             fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig("ks_kolmogorov_convergence.png", dpi=150, bbox_inches="tight")
plt.show()

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

np.random.seed(42)

# --- 2つの標本を生成 ---
n1, n2 = 50, 60
sample1 = np.random.normal(0, 1, n1)         # N(0, 1)
sample2 = np.random.normal(0.5, 1.2, n2)     # N(0.5, 1.2)

# SciPyでの二標本KS検定
D_stat, p_val = stats.ks_2samp(sample1, sample2)

# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))

# EDF の比較
ax = axes[0]
x_all = np.sort(np.concatenate([sample1, sample2]))
x_plot = np.concatenate([[-4], x_all, [5]])

# 各標本のEDF
edf1 = np.array([np.mean(sample1 <= x) for x in x_plot])
edf2 = np.array([np.mean(sample2 <= x) for x in x_plot])

ax.step(x_plot, edf1, where="post", linewidth=2, color="blue",
        label=f"Sample 1 (n={n1})")
ax.step(x_plot, edf2, where="post", linewidth=2, color="orange",
        label=f"Sample 2 (n={n2})")

# 最大乖離の表示
diff = np.abs(edf1 - edf2)
max_idx = np.argmax(diff)
x_max = x_plot[max_idx]
ax.annotate("", xy=(x_max, edf1[max_idx]),
            xytext=(x_max, edf2[max_idx]),
            arrowprops=dict(arrowstyle="<->", color="green", lw=2))
ax.text(x_max + 0.15, (edf1[max_idx] + edf2[max_idx]) / 2,
        f"D = {D_stat:.3f}", fontsize=11, color="green", fontweight="bold")

ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("Cumulative probability", fontsize=12)
ax.set_title(f"Two-sample KS test (p = {p_val:.4f})", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# ヒストグラム
ax = axes[1]
bins = np.linspace(-4, 5, 30)
ax.hist(sample1, bins=bins, density=True, alpha=0.5, color="blue",
        edgecolor="white", label=f"Sample 1: N(0, 1)")
ax.hist(sample2, bins=bins, density=True, alpha=0.5, color="orange",
        edgecolor="white", label=f"Sample 2: N(0.5, 1.2)")
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("Density", fontsize=12)
ax.set_title("Data distributions", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("ks_test_two_sample.png", dpi=150, bbox_inches="tight")
plt.show()

print(f"二標本KS統計量: D = {D_stat:.4f}")
print(f"p値: {p_val:.6f}")

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

np.random.seed(42)

n = 50
n_sim = 3000
alpha = 0.05

# 3つのタイプの対立仮説
# (1) 位置の差、(2) スケールの差、(3) 形状の差
shift_vals = np.linspace(0, 1.5, 12)

power_ks_shift = []
power_t_shift = []
power_ks_scale = []
power_ks_shape = []

for delta in shift_vals:
    rej_ks_s, rej_t_s, rej_ks_sc, rej_ks_sh = 0, 0, 0, 0

    for _ in range(n_sim):
        x_base = np.random.randn(n)

        # (1) 位置シフト
        y_shift = np.random.randn(n) + delta
        _, p_ks = stats.ks_2samp(x_base, y_shift)
        _, p_t = stats.ttest_ind(x_base, y_shift)
        if p_ks < alpha: rej_ks_s += 1
        if p_t < alpha: rej_t_s += 1

        # (2) スケール変化
        y_scale = np.random.randn(n) * (1 + delta)
        _, p_ks = stats.ks_2samp(x_base, y_scale)
        if p_ks < alpha: rej_ks_sc += 1

        # (3) 形状変化（混合）
        y_shape = np.where(
            np.random.rand(n) < 0.5 + delta / 3,
            np.random.randn(n),
            np.random.randn(n) + 3
        )
        _, p_ks = stats.ks_2samp(x_base, y_shape)
        if p_ks < alpha: rej_ks_sh += 1

    power_ks_shift.append(rej_ks_s / n_sim)
    power_t_shift.append(rej_t_s / n_sim)
    power_ks_scale.append(rej_ks_sc / n_sim)
    power_ks_shape.append(rej_ks_sh / n_sim)

fig, axes = plt.subplots(1, 3, figsize=(16, 5))

ax = axes[0]
ax.plot(shift_vals, power_ks_shift, "b-o", markersize=5, label="KS test")
ax.plot(shift_vals, power_t_shift, "r-s", markersize=5, label="t-test")
ax.axhline(alpha, color="gray", linestyle="--", linewidth=0.8)
ax.set_title("Location shift", fontsize=13)
ax.set_xlabel("Shift amount", fontsize=11)
ax.set_ylabel("Power", fontsize=11)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 1.05)

ax = axes[1]
ax.plot(shift_vals, power_ks_scale, "b-o", markersize=5, label="KS test")
ax.axhline(alpha, color="gray", linestyle="--", linewidth=0.8)
ax.set_title("Scale change", fontsize=13)
ax.set_xlabel("Scale increase", fontsize=11)
ax.set_ylabel("Power", fontsize=11)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 1.05)

ax = axes[2]
ax.plot(shift_vals, power_ks_shape, "b-o", markersize=5, label="KS test")
ax.axhline(alpha, color="gray", linestyle="--", linewidth=0.8)
ax.set_title("Shape change (mixture)", fontsize=13)
ax.set_xlabel("Mixture proportion change", fontsize=11)
ax.set_ylabel("Power", fontsize=11)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 1.05)

plt.tight_layout()
plt.savefig("ks_test_power.png", dpi=150, bbox_inches="tight")
plt.show()