Q-Q 그림 신뢰띠 모의실험¶

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

def qq_with_band(x, B=800, seed=42):
    """Q-Q 그림에 모의실험으로 만든 95% 띠를 얹는다.

    Q-Q 그림의 점들은 웬만큼 흔들리기 마련이라, 직선에서 조금 벗어난 것이
    문제인지 아닌지 눈으로는 알기 어렵다. 적합한 정규분포에서 같은 크기의
    표본을 B 번 뽑아 각 자리의 2.5·97.5 백분위점을 구하면, "정규라면 이
    정도까지는 흔들린다"는 범위를 그릴 수 있다.
    """
    x = np.asarray(x, dtype=float)
    n = x.size
    mu, sd = x.mean(), x.std(ddof=1)

    # 이론 분위수. 0.5 를 빼는 것은 i/n 이 마지막 점에서 1 이 되어
    # ppf 가 무한이 되는 것을 피하기 위한 흔한 보정이다.
    p = (np.arange(1, n + 1) - 0.5) / n
    q_theor = stats.norm.ppf(p)

    x_sorted = np.sort(x)

    # 적합된 정규분포에서 B 번 표본을 뽑아 각 순서통계량의 분포를 얻는다.
    rng = np.random.default_rng(seed)
    sims = np.sort(rng.normal(mu, sd, size=(B, n)), axis=1)
    lo = np.percentile(sims, 2.5, axis=0)
    hi = np.percentile(sims, 97.5, axis=0)

    fig, ax = plt.subplots(figsize=(7, 4))
    ax.scatter(q_theor, x_sorted, s=15)
    ax.plot(q_theor, mu + sd * q_theor, linestyle="--")
    ax.fill_between(q_theor, lo, hi, alpha=0.15,
                    label="95% pointwise band")
    ax.set_title("Q-Q Plot with Simulated 95% Band")
    ax.set_xlabel("Theoretical quantiles (Normal)")
    ax.set_ylabel("Ordered data")
    ax.legend()
    plt.tight_layout()
    plt.show()

rng = np.random.default_rng(123)
x = rng.lognormal(mean=0.0, sigma=0.6, size=300)
qq_with_band(x, B=600, seed=7)

def runs(mask):
    """True 인 자리를 연속 구간으로 묶는다 (1 부터 센 순위로)."""
    idx = np.where(mask)[0]
    if idx.size == 0:
        return []
    out, start, prev = [], idx[0], idx[0]
    for j in idx[1:]:
        if j == prev + 1:
            prev = j
        else:
            out.append((start + 1, prev + 1))
            start = prev = j
    out.append((start + 1, prev + 1))
    return out

rng = np.random.default_rng(123)
x = rng.lognormal(mean=0.0, sigma=0.6, size=300)

n = x.size
mu, sd = x.mean(), x.std(ddof=1)
p = (np.arange(1, n + 1) - 0.5) / n
q = stats.norm.ppf(p)
x_sorted = np.sort(x)

sims = np.sort(np.random.default_rng(7).normal(mu, sd, size=(600, n)), axis=1)
lo = np.percentile(sims, 2.5, axis=0)
hi = np.percentile(sims, 97.5, axis=0)

above, below = x_sorted > hi, x_sorted < lo
print(f"적합된 정규: mu-hat {mu:.4f}, sd-hat {sd:.4f}")
print(f"띠 밖의 점 {above.sum() + below.sum()}/{n} "
      f"({100 * (above.sum() + below.sum()) / n:.1f}%)")
print(f"  위로 벗어난 구간(순위): {runs(above)}")
print(f"  아래로 벗어난 순위 범위: {np.where(below)[0].min() + 1}"
      f" ~ {np.where(below)[0].max() + 1}  ({below.sum()}점)")

width = hi - lo
print(f"띠 폭: 가운데 최솟값 {width.min():.4f}, 양 끝 {width[0]:.4f} / "
      f"{width[-1]:.4f}, 최대/최소 = {width.max() / width.min():.2f}")
print(f"띠 하한의 최솟값 {lo.min():.4f}   자료 최솟값 {x_sorted[0]:.4f}")
print(f"P(적합된 정규 < 0) = {stats.norm.cdf(0, mu, sd):.4f}")
print(f"최댓값 {x_sorted[-1]:.4f}  띠 상한 {hi[-1]:.4f}  "
      f"({x_sorted[-1] / hi[-1]:.2f}배)")

s2 = 0.6 ** 2
print(f"\n이론: 평균 {np.exp(s2 / 2):.4f}  sd "
      f"{np.sqrt((np.exp(s2) - 1) * np.exp(s2)):.4f}  왜도 "
      f"{(np.exp(s2) + 2) * np.sqrt(np.exp(s2) - 1):.4f}")
print(f"표본: 평균 {mu:.4f}  sd {sd:.4f}  g1 {stats.skew(x):.4f}"
      f"  G1 {stats.skew(x, bias=False):.4f}")
print(f"샤피로-윌크 p = {stats.shapiro(x)[1]:.4g}")

적합된 정규: mu-hat 1.2044, sd-hat 0.7138
띠 밖의 점 204/300 (68.0%)
  위로 벗어난 구간(순위): [(1, 48), (278, 293), (295, 300)]
  아래로 벗어난 순위 범위: 103 ~ 241  (134점)
띠 폭: 가운데 최솟값 0.1994, 양 끝 1.0862 / 1.0342, 최대/최소 = 5.45
띠 하한의 최솟값 -1.4917   자료 최솟값 0.1807
P(적합된 정규 < 0) = 0.0458
최댓값 4.8422  띠 상한 3.8469  (1.26배)

이론: 평균 1.1972  sd 0.7881  왜도 2.2601
표본: 평균 1.2044  sd 0.7138  g1 1.5056  G1 1.5132
샤피로-윌크 p = 3.877e-14

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

rng = np.random.default_rng(0)
x = rng.normal(0, 1, size=200)

n = x.size
mu, sd = x.mean(), x.std(ddof=1)
p = (np.arange(1, n + 1) - 0.5) / n
q = stats.norm.ppf(p)
x_sorted = np.sort(x)

sims = np.sort(rng.normal(mu, sd, size=(1000, n)), axis=1)
lo = np.percentile(sims, 2.5, axis=0)
hi = np.percentile(sims, 97.5, axis=0)

outside = np.sum((x_sorted < lo) | (x_sorted > hi))
print(f"Points outside band: {outside} / {n}")

fig, ax = plt.subplots(figsize=(7, 4))
ax.scatter(q, x_sorted, s=15)
ax.plot(q, mu + sd * q, linestyle="--")
ax.fill_between(q, lo, hi, alpha=0.15, label="95% band")
ax.legend()
ax.set_title("Q-Q Plot with 95% Band (Normal Data)")
plt.tight_layout()
plt.show()

Points outside band: 1 / 200

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

rng = np.random.default_rng(1)
x = rng.standard_t(df=4, size=200)

n = x.size
mu, sd = x.mean(), x.std(ddof=1)
p = (np.arange(1, n + 1) - 0.5) / n
q = stats.norm.ppf(p)
x_sorted = np.sort(x)

sims = np.sort(rng.normal(mu, sd, size=(1000, n)), axis=1)
lo = np.percentile(sims, 2.5, axis=0)
hi = np.percentile(sims, 97.5, axis=0)

fig, ax = plt.subplots(figsize=(7, 4))
ax.scatter(q, x_sorted, s=15)
ax.plot(q, mu + sd * q, linestyle="--")
ax.fill_between(q, lo, hi, alpha=0.15, label="95% band")
ax.legend()
ax.set_title("Q-Q Plot with 95% Band (t4 Data)")
plt.tight_layout()
plt.show()

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

rng = np.random.default_rng(0)
x = rng.normal(0, 1, size=200)

n = x.size
mu, sd = x.mean(), x.std(ddof=1)
p = (np.arange(1, n + 1) - 0.5) / n
q = stats.norm.ppf(p)
x_sorted = np.sort(x)

sims = np.sort(rng.normal(mu, sd, size=(1000, n)), axis=1)
lo95 = np.percentile(sims, 2.5, axis=0)
hi95 = np.percentile(sims, 97.5, axis=0)
lo99 = np.percentile(sims, 0.5, axis=0)
hi99 = np.percentile(sims, 99.5, axis=0)

print(f"Median width ratio: {np.median((hi99 - lo99) / (hi95 - lo95)):.3f}")

fig, ax = plt.subplots(figsize=(7, 4))
ax.fill_between(q, lo99, hi99, alpha=0.10, label="99% band")
ax.fill_between(q, lo95, hi95, alpha=0.15, label="95% band")
ax.scatter(q, x_sorted, s=15, zorder=3)
ax.plot(q, mu + sd * q, linestyle="--")
ax.legend()
ax.set_title("95% vs 99% Pointwise Bands")
plt.tight_layout()
plt.show()

Median width ratio: 1.305