상관 검정 개관¶

['ipykernel_launcher.py', '-f', '/tmp/kernel-abc.json']

import argparse
import numpy as np

parser = argparse.ArgumentParser(description='Correlation Test Examples')
parser.add_argument('--seed', type=int, default=1, metavar='S',
                    help='random seed (default: 1)')
# parse_args()가 아니라 parse_known_args()를 쓴다.
# 노트북이나 REPL에서는 sys.argv에 다른 인자가 들어 있어 parse_args()가
# SystemExit을 던지며 커널을 멈춰 세운다.
ARGS, _ = parser.parse_known_args()

np.random.seed(ARGS.seed)
ARGS.size = 1000

# (1) 노트북이 남기는 찌꺼기 인자를 흉내 내어 두 방식을 견준다.
import sys
import contextlib
import io

saved = sys.argv
sys.argv = ['ipykernel_launcher.py', '-f', '/tmp/kernel-abc.json']
try:
    a, unknown = parser.parse_known_args()
    print(f"parse_known_args: 성공, seed = {a.seed}, 버린 인자 = {unknown}")
    with contextlib.redirect_stderr(io.StringIO()) as err:
        try:
            parser.parse_args()
            print("parse_args: 성공")
        except SystemExit as ex:
            msg = err.getvalue().strip().splitlines()[-1]
            print(f"parse_args: SystemExit(code={ex.code}) — "
                  f"error: {msg.split('error: ', 1)[-1]}")
finally:
    sys.argv = saved

# (2) seed 가 고정하는 것과 고정하지 않는 것
np.random.seed(1)
a1 = np.random.rand(3)
np.random.seed(1)
a2 = np.random.rand(3)
print(f"\nnp.random.rand  1회차 {a1.round(6)}")
print(f"np.random.rand  2회차 {a2.round(6)}   같은가: {np.array_equal(a1, a2)}")

np.random.seed(1)
b1 = np.random.default_rng().random(3)
np.random.seed(1)
b2 = np.random.default_rng().random(3)
print(f"default_rng()   두 번 뽑아 같은가: {np.array_equal(b1, b2)}   "
      f"(돌릴 때마다 값이 달라지므로 값은 적지 않는다)")

np.random.seed(1)
c1 = np.random.default_rng(1).random(3)
c2 = np.random.default_rng(1).random(3)
print(f"default_rng(1)  두 번  {c1.round(6)}   같은가: {np.array_equal(c1, c2)}")
print(f"  (np.random.seed(1) 을 걸어도 default_rng(1) 의 값은 그것과 무관하다)")

parse_known_args: 성공, seed = 1, 버린 인자 = ['-f', '/tmp/kernel-abc.json']
parse_args: SystemExit(code=2) — error: unrecognized arguments: -f /tmp/kernel-abc.json

np.random.rand  1회차 [4.17022e-01 7.20324e-01 1.14000e-04]
np.random.rand  2회차 [4.17022e-01 7.20324e-01 1.14000e-04]   같은가: True
default_rng()   두 번 뽑아 같은가: False   (돌릴 때마다 값이 달라지므로 값은 적지 않는다)
default_rng(1)  두 번  [0.511822 0.950464 0.14416 ]   같은가: True
  (np.random.seed(1) 을 걸어도 default_rng(1) 의 값은 그것과 무관하다)

import numpy as np
from scipy import stats

# 위의 global_name_space.py를 파일로 저장했다면 다음 한 줄로 대신할 수 있다.
#   from global_name_space import ARGS

def load_data(data_type=0):
    data_dict = {}

    x = np.random.rand(ARGS.size) * 20
    eps = np.random.rand(ARGS.size) * 10

    # 자료 0: 관계 없음
    y = np.random.rand(ARGS.size) * 20
    data_dict[0] = (x, y)

    # 자료 1: 단조이지만 곡선인 관계(삼차)
    y = (x + eps) ** 3
    data_dict[1] = (x, y)

    # 자료 2: 단조가 아닌 관계(사인)
    y = np.sin(x + eps)
    data_dict[2] = (x, y)

    return data_dict

np.random.seed(1)
class ARGS: size = 1000
d = load_data()

# (1) 세 자료가 x 를 공유하는가
print(f"자료 0 과 1 의 x 가 같은 배열인가: {np.array_equal(d[0][0], d[1][0])}")
print(f"자료 1 과 2 의 x 가 같은 배열인가: {np.array_equal(d[1][0], d[2][0])}")

# (2) 자료 1 의 y 가 x 의 단조함수인가 — 반례를 찾는다
x, y = d[1]
o = np.argsort(x)
xs, ys = x[o], y[o]
bad = np.flatnonzero(np.diff(ys) < 0)
print(f"\n자료 1: x 를 오름차순으로 놓았을 때 y 가 줄어드는 자리 = "
      f"{len(bad)} 곳 / {len(x)-1}")
i = bad[0]
print(f"  예: x = {xs[i]:.4f} -> {xs[i+1]:.4f} 인데 y = {ys[i]:.1f} -> {ys[i+1]:.1f}")
print(f"  (x 가 커져도 eps 가 더 작으면 x+eps 가 줄 수 있다)")

# 세제곱은 순위를 바꾸지 않는다. y = (x+eps)^3 이므로 cbrt(y) = x+eps 다.
print(f"\n세제곱은 순위를 바꾸지 않는다:")
print(f"  spearman(x, (x+eps)^3) = {stats.spearmanr(x, y).statistic:.6f}")
print(f"  spearman(x,  x+eps   ) = {stats.spearmanr(x, np.cbrt(y)).statistic:.6f}")
print(f"  pearson (x, (x+eps)^3) = {stats.pearsonr(x, y).statistic:.6f}")
print(f"  pearson (x,  x+eps   ) = {stats.pearsonr(x, np.cbrt(y)).statistic:.6f}")

# (3) 자료 2 가 도는 주기 수
u2 = np.cbrt(d[1][1])          # = x + eps
print(f"\n자료 2: x+eps 의 범위 = [{u2.min():.4f}, {u2.max():.4f}],  "
      f"sin 이 도는 주기 = {(u2.max() - u2.min()) / (2 * np.pi):.3f}")

자료 0 과 1 의 x 가 같은 배열인가: True
자료 1 과 2 의 x 가 같은 배열인가: True

자료 1: x 를 오름차순으로 놓았을 때 y 가 줄어드는 자리 = 503 곳 / 999
  예: x = 0.0023 -> 0.0080 인데 y = 425.2 -> 21.1
  (x 가 커져도 eps 가 더 작으면 x+eps 가 줄 수 있다)

세제곱은 순위를 바꾸지 않는다:
  spearman(x, (x+eps)^3) = 0.900004
  spearman(x,  x+eps   ) = 0.900004
  pearson (x, (x+eps)^3) = 0.815408
  pearson (x,  x+eps   ) = 0.894142

자료 2: x+eps 의 범위 = [0.5744, 29.6839],  sin 이 도는 주기 = 4.633

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

# 위의 load_data.py를 파일로 저장했다면 다음 한 줄로 대신할 수 있다.
#   from load_data import load_data

def main():
    data_dict = load_data()
    _, axes = plt.subplots(1, len(data_dict), figsize=(12, 3))

    for ax, (x, y) in zip(axes, data_dict.values()):
        ax.plot(x, y, ".k")
        coef, p_val = stats.pearsonr(x, y)
        ax.set_title(f"Pearson's r: {coef:.4f}\np-value: {p_val:.4f}")
    plt.tight_layout()
    plt.show()

np.random.seed(1)          # 보기 1 의 설정 모듈이 하는 일 (별도 실행이면 자동으로 된다)

if __name__ == "__main__":
    main()

# 모집단 값을 적률로 정확히 구한다.
from fractions import Fraction as F
from math import comb

# E[X^k] = 20^k/(k+1),  E[eps^k] = 10^k/(k+1)
EX = [F(20**k, k + 1) for k in range(7)]
EE = [F(10**k, k + 1) for k in range(7)]
EU = lambda m: sum(comb(m, k) * EX[m - k] * EE[k] for k in range(m + 1))

EXU3 = EX[4] + 3 * EX[3] * EE[1] + 3 * EX[2] * EE[2] + EX[1] * EE[3]
cov = EXU3 - EX[1] * EU(3)
varX = EX[2] - EX[1] ** 2
varU3 = EU(6) - EU(3) ** 2
rho1 = float(cov) / np.sqrt(float(varX) * float(varU3))
print("자료 1 의 모집단 Pearson 상관")
print(f"  E[X U^3] = {EXU3} = {float(EXU3):.4f}")
print(f"  E[U^3]   = {EU(3)} = {float(EU(3)):.4f}")
print(f"  Cov      = {cov} = {float(cov):.4f}")
print(f"  Var(X)   = {varX} = {float(varX):.4f}")
print(f"  E[U^6]   = {EU(6)} = {float(EU(6)):.4f}")
print(f"  Var(U^3) = {varU3} = {float(varU3):.4f}")
print(f"  rho      = {rho1:.6f}")

# 자료 2 는 사인이 들어가 적률로 떨어지지 않으므로 수치적분으로 구한다.
from scipy import integrate
Es = integrate.dblquad(lambda e, x: np.sin(x + e) / 200, 0, 20, 0, 10)[0]
Exs = integrate.dblquad(lambda e, x: x * np.sin(x + e) / 200, 0, 20, 0, 10)[0]
Es2 = integrate.dblquad(lambda e, x: np.sin(x + e) ** 2 / 200, 0, 20, 0, 10)[0]
rho2 = (Exs - 10 * Es) / np.sqrt(float(varX) * (Es2 - Es ** 2))
print(f"\n자료 2 의 모집단 Pearson 상관 = {rho2:.6f}   (정확히 0 이 아니다)")

np.random.seed(1)          # 그림이 쓴 것과 같은 자료를 다시 얻는다
d = load_data()
print(f"\n{'자료':>4s} {'표본 r':>10s} {'모집단 rho':>11s} {'SE':>8s} {'z':>7s} {'p':>8s}")
pop = [0.0, rho1, rho2]
for k, (x, y) in d.items():
    r, p = stats.pearsonr(x, y)
    se = (1 - pop[k] ** 2) / np.sqrt(ARGS.size)
    print(f"{k:4d} {r:+10.6f} {pop[k]:+11.6f} {se:8.6f} {(r - pop[k]) / se:+7.3f} {p:8.4f}")

자료 1 의 모집단 Pearson 상관
  E[X U^3] = 233500/3 = 77833.3333
  E[U^3]   = 5250 = 5250.0000
  Cov      = 76000/3 = 25333.3333
  Var(X)   = 100/3 = 33.3333
  E[U^6]   = 394000000/7 = 56285714.2857
  Var(U^3) = 201062500/7 = 28723214.2857
  rho      = 0.818722

자료 2 의 모집단 Pearson 상관 = 0.027999   (정확히 0 이 아니다)

  자료       표본 r     모집단 rho       SE       z        p
   0  +0.034844   +0.000000 0.031623  +1.102   0.2710
   1  +0.815408   +0.818722 0.010426  -0.318   0.0000
   2  -0.018298   +0.027999 0.031598  -1.465   0.5633

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

# 위의 load_data.py를 파일로 저장했다면 다음 한 줄로 대신할 수 있다.
#   from load_data import load_data

def main():
    data_dict = load_data()
    _, axes = plt.subplots(1, len(data_dict), figsize=(12, 3))

    for ax, (x, y) in zip(axes, data_dict.values()):
        ax.plot(x, y, ".k")
        coef, p_val = stats.spearmanr(x, y)
        ax.set_title(f"Spearman rho: {coef:.4f}\np-value: {p_val:.4f}")
    plt.tight_layout()
    plt.show()

np.random.seed(1)          # 보기 1 의 설정 모듈이 하는 일 (별도 실행이면 자동으로 된다)

if __name__ == "__main__":
    main()

# 복사 항등식을 기호적분으로 확인한다.
import sympy as sp

x_, e_ = sp.symbols('x e', real=True)
# U = X + eps 의 분포함수 (U(0,20) 과 U(0,10) 의 합성곱, 사다리꼴)
u_ = x_ + e_
F_U = sp.Piecewise((u_**2 / 400, u_ <= 10),
                   (sp.Rational(1, 4) + (u_ - 10) / 20, u_ <= 20),
                   (1 - (30 - u_)**2 / 400, True))
# E[F_X(X) F_U(U)],  F_X(x) = x/20,  결합밀도 1/200
EFF = sp.integrate(sp.integrate((x_ / 20) * F_U / 200, (e_, 0, 10)), (x_, 0, 20))
EFF = sp.nsimplify(sp.simplify(EFF))
rho_s = sp.simplify(12 * EFF - 3)
print(f"E[F_X(X) F_U(U)] = {EFF} = {float(EFF):.6f}")
print(f"rho_s = 12 E[...] - 3 = {rho_s} = {float(rho_s):.6f}")

# 아주 큰 표본으로 확인
rng = np.random.default_rng(99)
vals = [stats.spearmanr(a := rng.random(200_000) * 20,
                        (a + rng.random(200_000) * 10) ** 3).statistic
        for _ in range(20)]
print(f"모의 rho_s = {np.mean(vals):.6f} +- {np.std(vals) / np.sqrt(20):.6f}")

# 자료 2 의 모집단 rho_s 는 닫힌 꼴이 없으므로 모의로 잰다
vals2 = [stats.spearmanr(a := rng.random(200_000) * 20,
                         np.sin(a + rng.random(200_000) * 10)).statistic
         for _ in range(20)]
rho_s2 = np.mean(vals2)
print(f"자료 2 모의 rho_s = {rho_s2:.6f} +- {np.std(vals2) / np.sqrt(20):.6f}")

np.random.seed(1)          # 그림이 쓴 것과 같은 자료를 다시 얻는다
d = load_data()
pop = [0.0, float(rho_s), rho_s2]
print(f"\n{'자료':>4s} {'표본 r_s':>10s} {'모집단':>10s} {'SE':>8s} {'z':>7s}")
for k, (x, y) in d.items():
    rs = stats.spearmanr(x, y).statistic
    se = (1 - pop[k] ** 2) / np.sqrt(ARGS.size)
    print(f"{k:4d} {rs:+10.6f} {pop[k]:+10.6f} {se:8.6f} {(rs - pop[k]) / se:+7.3f}")

E[F_X(X) F_U(U)] = 13/40 = 0.325000
rho_s = 12 E[...] - 3 = 9/10 = 0.900000
모의 rho_s = 0.900015 +- 0.000088
자료 2 모의 rho_s = 0.027899 +- 0.000549

  자료     표본 r_s        모집단       SE       z
   0  +0.036465  +0.000000 0.031623  +1.153
   1  +0.900004  +0.900000 0.006008  +0.001
   2  -0.017612  +0.027899 0.031598  -1.440

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
from fractions import Fraction as F

# 위의 load_data.py를 파일로 저장했다면 다음 한 줄로 대신할 수 있다.
#   from load_data import load_data

def main():
    data_dict = load_data()
    _, axes = plt.subplots(1, len(data_dict), figsize=(12, 3))

    for ax, (x, y) in zip(axes, data_dict.values()):
        ax.plot(x, y, ".k")
        coef, p_val = stats.kendalltau(x, y)
        ax.set_title(f"Kendall's tau: {coef:.4f}\np-value: {p_val:.4f}")
    plt.tight_layout()
    plt.show()

np.random.seed(1)          # 보기 1 의 설정 모듈이 하는 일 (별도 실행이면 자동으로 된다)

if __name__ == "__main__":
    main()

tau_exact = F(17, 24)
print(f"자료 1 의 모집단 tau = {tau_exact} = {float(tau_exact):.6f}")

# 유도의 중간값들을 수치적분으로 확인한다.
from scipy import integrate
fD = lambda d: (20 - abs(d)) / 400          # D = X1 - X2 의 밀도 (삼각)
fG = lambda g: (10 - abs(g)) / 100          # G = e1 - e2 의 밀도 (삼각)
PG = lambda d: integrate.quad(fG, -10, -d)[0]
half = integrate.quad(lambda d: fD(d) * PG(d), 0, 10)[0]
print(f"  P(D>0, D+G<0) = {half:.8f}   손계산 7/96 = {float(F(7,96)):.8f}")
print(f"  P(부호 뒤집힘) = {2*half:.8f}   손계산 7/48 = {float(F(7,48)):.8f}")
print(f"  tau = 1 - 2*{2*half:.6f} = {1 - 4*half:.6f}")

rng = np.random.default_rng(99)
vals = [stats.kendalltau(a := rng.random(100_000) * 20,
                         (a + rng.random(100_000) * 10) ** 3).statistic
        for _ in range(20)]
print(f"  모의 tau = {np.mean(vals):.6f} +- {np.std(vals) / np.sqrt(20):.6f}")

vals2 = [stats.kendalltau(a := rng.random(100_000) * 20,
                          np.sin(a + rng.random(100_000) * 10)).statistic
         for _ in range(20)]
tau2 = np.mean(vals2)
print(f"\n자료 2 모의 tau = {tau2:.6f} +- {np.std(vals2) / np.sqrt(20):.6f}")

n = ARGS.size
se0 = np.sqrt(2 * (2 * n + 5) / (9 * n * (n - 1)))   # H0 아래 tau 의 표준편차
print(f"\nH0 아래 sd(tau) = sqrt(2(2n+5)/(9n(n-1))) = {se0:.6f}")

np.random.seed(1)          # 그림이 쓴 것과 같은 자료를 다시 얻는다
d = load_data()
pop = [0.0, float(tau_exact), tau2]
print(f"\n{'자료':>4s} {'표본 tau':>10s} {'모집단':>10s} {'비교용 rho_s':>12s}")
for k, (x, y) in d.items():
    t = stats.kendalltau(x, y).statistic
    rs = stats.spearmanr(x, y).statistic
    print(f"{k:4d} {t:+10.6f} {pop[k]:+10.6f} {rs:+12.6f}")

자료 1 의 모집단 tau = 17/24 = 0.708333
  P(D>0, D+G<0) = 0.07291667   손계산 7/96 = 0.07291667
  P(부호 뒤집힘) = 0.14583333   손계산 7/48 = 0.14583333
  tau = 1 - 2*0.145833 = 0.708333
  모의 tau = 0.708444 +- 0.000187

자료 2 모의 tau = 0.018885 +- 0.000400

H0 아래 sd(tau) = sqrt(2(2n+5)/(9n(n-1))) = 0.021119

  자료     표본 tau        모집단    비교용 rho_s
   0  +0.023764  +0.000000    +0.036465
   1  +0.709373  +0.708333    +0.900004
   2  -0.012344  +0.018885    -0.017612

나이 : 18  24  25  26  33  37  40  45  57  64
소득 : 15  26  29  32  33  41  45  52  68  80   (천 달러)

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

def main():
    x = [18, 25, 57, 45, 26, 64, 37, 40, 24, 33]
    y = [15_000, 29_000, 68_000, 52_000, 32_000, 80_000, 41_000, 45_000, 26_000, 33_000]

    coef, p_val = stats.pearsonr(x, y)
    print(f"Pearson's r:   coef = {coef:.4f},  p-value = {p_val:.4f}")

    coef, p_val = stats.spearmanr(x, y)
    print(f"Spearman rho:  coef = {coef:.4f},  p-value = {p_val:.4f}")

    coef, p_val = stats.kendalltau(x, y)
    print(f"Kendall's tau: coef = {coef:.4f},  p-value = {p_val:.4f}")

    fig, ax = plt.subplots(figsize=(12, 3))
    ax.plot(x, y, "ok")
    ax.set_xlabel("Age")
    ax.set_ylabel("Income")
    plt.show()

if __name__ == "__main__":
    main()

x = np.array([18, 25, 57, 45, 26, 64, 37, 40, 24, 33])
y = np.array([15_000, 29_000, 68_000, 52_000, 32_000, 80_000,
              41_000, 45_000, 26_000, 33_000])
n = len(x)

# (1) 순위가 완전히 맞아떨어지는가
o = np.argsort(x)
print(f"\n나이 오름차순 : {x[o].tolist()}")
print(f"그때의 소득    : {y[o].tolist()}")
print(f"소득도 오름차순인가: {bool(np.all(np.diff(y[o]) > 0))}")

# 지수 관계인가 — 세 가지 눈금에서 Pearson 을 재 본다
print(f"\nr(x, y)         = {stats.pearsonr(x, y).statistic:.6f}   (직선)")
print(f"r(x, log y)     = {stats.pearsonr(x, np.log(y)).statistic:.6f}   (지수)")
print(f"r(log x, log y) = {stats.pearsonr(np.log(x), np.log(y)).statistic:.6f}   (거듭제곱)")

# (2) 세 p-값
print(f"\nPearson  p = {stats.pearsonr(x, y).pvalue:.6e}")
rs = stats.spearmanr(x, y)
print(f"Spearman p = {rs.pvalue:.6e}   <- 이 값은 무엇인가")
print(f"  r_s = {rs.statistic!r}")
print(f"  1 - r_s^2 = {1 - rs.statistic**2:.6e}  (0 이어야 하는데 반올림 오차가 남는다)")
print(f"  t = r_s sqrt((n-2)/(1-r_s^2)) = "
      f"{rs.statistic * np.sqrt((n-2)/(1-rs.statistic**2)):.6e}")
kt = stats.kendalltau(x, y)
print(f"Kendall  p = {kt.pvalue:.6e}")
exact = 2 / factorial(n)
print(f"\n순위가 완전히 맞을 확률 (정확 순열): 2/{n}! = {exact:.6e}")
print(f"  Kendall 의 p 와 같은가: {kt.pvalue == exact}")
print(f"  Spearman 의 p 와의 비: {rs.pvalue / exact:.3e}")

Pearson's r:   coef = 0.9923,  p-value = 0.0000
Spearman rho:  coef = 1.0000,  p-value = 0.0000
Kendall's tau: coef = 1.0000,  p-value = 0.0000

나이 오름차순 : [18, 24, 25, 26, 33, 37, 40, 45, 57, 64]
그때의 소득    : [15000, 26000, 29000, 32000, 33000, 41000, 45000, 52000, 68000, 80000]
소득도 오름차순인가: True

r(x, y)         = 0.992285   (직선)
r(x, log y)     = 0.957012   (지수)
r(log x, log y) = 0.982551   (거듭제곱)

Pearson  p = 1.535456e-08
Spearman p = 6.646897e-64   <- 이 값은 무엇인가
  r_s = 0.9999999999999999
  1 - r_s^2 = 2.220446e-16  (0 이어야 하는데 반올림 오차가 남는다)
  t = r_s sqrt((n-2)/(1-r_s^2)) = 1.898125e+08
Kendall  p = 5.511464e-07

순위가 완전히 맞을 확률 (정확 순열): 2/10! = 5.511464e-07
  Kendall 의 p 와 같은가: True
  Spearman 의 p 와의 비: 1.206e-57

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(30001)
B = 20_000
gens = {
    "정규": lambda n: rng.standard_normal(n),
    "균등": lambda n: rng.uniform(-1, 1, n),
    "지수 (왜도 2)": lambda n: rng.exponential(1, n),
    "로그정규 (왜도 6)": lambda n: np.exp(rng.standard_normal(n)),
    "t(3) 두꺼운 꼬리": lambda n: rng.standard_t(3, n),
    "코시": lambda n: rng.standard_cauchy(n),
}
print("두 변수가 독립일 때의 1종 오류 (명목 0.05)")
print(f"{'분포':>22s} {'n=10':>8s} {'n=30':>8s} {'n=100':>8s}")
for lab, g in gens.items():
    row = []
    for n in [10, 30, 100]:
        a = sum(stats.pearsonr(g(n), g(n)).pvalue < 0.05 for _ in range(B))
        row.append(a / B)
    print(f"{lab:>22s} {row[0]:8.4f} {row[1]:8.4f} {row[2]:8.4f}")

두 변수가 독립일 때의 1종 오류 (명목 0.05)
                    분포     n=10     n=30    n=100
                    정규   0.0498   0.0483   0.0500
                    균등   0.0498   0.0500   0.0534
             지수 (왜도 2)   0.0483   0.0530   0.0522
           로그정규 (왜도 6)   0.0578   0.0530   0.0451
           t(3) 두꺼운 꼬리   0.0474   0.0525   0.0505
                    코시   0.0771   0.0641   0.0470

독립 + H0 → (x_i, y_σ(i)) 의 분포가 모든 순열 σ 에서 동일
         → r 의 귀무분포가 주변분포에 (거의) 의존하지 않는다
         → t 근사가 잘 맞는다

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(30002)

def gen(kind, n, rho=0.6):
    """가우스 코퓰러로 주변분포만 바꾼다."""
    z = rng.standard_normal((n, 2))
    x = z[:, 0]
    y = rho * z[:, 0] + np.sqrt(1 - rho**2) * z[:, 1]
    if kind == "정규":
        return x, y
    if kind == "로그정규":
        return np.exp(x), np.exp(y)
    df = 5 if kind.startswith("t(5)") else 3
    return (stats.t.ppf(stats.norm.cdf(x), df),
            stats.t.ppf(stats.norm.cdf(y), df))

B = 10_000
print("참 ρ=0.6 의 피셔 z 구간, 명목 95%")
print(f"{'분포':>16s} {'n=30':>9s} {'n=100':>9s} {'n=500':>9s} {'참 r':>9s}")
for kind in ["정규", "로그정규", "t(5) 주변화", "t(3) 주변화"]:
    bx, by = gen(kind, 400_000)
    true = np.corrcoef(bx, by)[0, 1]
    row = []
    for n in [30, 100, 500]:
        a = 0
        for _ in range(B):
            x, y = gen(kind, n)
            r = np.corrcoef(x, y)[0, 1]
            q = 1.96 / np.sqrt(n - 3)
            a += (np.tanh(np.arctanh(r) - q) <= true
                  <= np.tanh(np.arctanh(r) + q))
        row.append(a / B)
    print(f"{kind:>16s} {row[0]:9.4f} {row[1]:9.4f} {row[2]:9.4f} {true:9.4f}")

참 ρ=0.6 의 피셔 z 구간, 명목 95%
              분포      n=30     n=100     n=500       참 r
              정규    0.9469    0.9486    0.9510    0.5988
            로그정규    0.7965    0.7332    0.6482    0.4765
        t(5) 주변화    0.9394    0.9380    0.9295    0.5862
        t(3) 주변화    0.9097    0.8697    0.7772    0.5372

rng = np.random.default_rng(30003)

def lognorm_pair(n, rho=0.6):
    z = rng.standard_normal((n, 2))
    return (np.exp(z[:, 0]),
            np.exp(rho * z[:, 0] + np.sqrt(1 - rho**2) * z[:, 1]))

bx, by = lognorm_pair(400_000)
true = np.corrcoef(bx, by)[0, 1]
print(f"\n로그정규, 참 상관 = {true:.4f}")
B, R = 2_000, 999
for n in [30, 100]:
    a = b = 0
    for _ in range(B):
        x, y = lognorm_pair(n)
        r = np.corrcoef(x, y)[0, 1]
        q = 1.96 / np.sqrt(n - 3)
        a += (np.tanh(np.arctanh(r) - q) <= true
              <= np.tanh(np.arctanh(r) + q))
        idx = rng.integers(0, n, (R, n))
        bs = np.array([np.corrcoef(x[i], y[i])[0, 1] for i in idx])
        lo, hi = np.quantile(bs, [0.025, 0.975])
        b += lo <= true <= hi
    print(f"  n={n:3d}: 피셔 z {a / B:.4f},  부트스트랩 백분위 {b / B:.4f}")

로그정규, 참 상관 = 0.4764
  n= 30: 피셔 z 0.7870,  부트스트랩 백분위 0.9110
  n=100: 피셔 z 0.7290,  부트스트랩 백분위 0.9135

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

def tost(r, n, delta=0.2):
    """상관의 두 단측 검정. 반환값이 작으면 등가."""
    z, se = np.arctanh(r), 1 / np.sqrt(n - 3)
    p1 = stats.norm.sf((z - np.arctanh(-delta)) / se)
    p2 = stats.norm.cdf((z - np.arctanh(delta)) / se)
    return max(p1, p2)

rng = np.random.default_rng(30003)
B = 10_000
print("등가 한계 ±0.2, 명목 0.05. 표는 '등가'라고 결론 내릴 확률")
print(f"{'참 ρ':>7s} {'n=50':>9s} {'n=200':>9s} {'n=500':>9s} {'n=1000':>9s}")
for rho in [0.0, 0.1, 0.2, 0.3]:
    row = []
    for n in [50, 200, 500, 1000]:
        a = 0
        for _ in range(B):
            z = rng.standard_normal((n, 2))
            r = np.corrcoef(z[:, 0], rho * z[:, 0]
                            + np.sqrt(1 - rho**2) * z[:, 1])[0, 1]
            a += tost(r, n) < 0.05
        row.append(a / B)
    print(f"{rho:7.1f} {row[0]:9.4f} {row[1]:9.4f} {row[2]:9.4f} {row[3]:9.4f}")

등가 한계 ±0.2, 명목 0.05. 표는 '등가'라고 결론 내릴 확률
    참 ρ      n=50     n=200     n=500    n=1000
    0.0    0.0000    0.7707    0.9966    1.0000
    0.1    0.0000    0.4097    0.7422    0.9411
    0.2    0.0000    0.0454    0.0507    0.0470
    0.3    0.0000    0.0007    0.0000    0.0000

90% 신뢰구간이 [-δ, +δ] 안에 완전히 들어가면 등가
  (양측 90% = 단측 5% 두 개)

이것이 TOST 와 정확히 같은 결정을 준다

r = 0.04,  95% CI [-0.10, 0.18],  n = 200
등가 검정 (δ = 0.2): p = 0.021  → 등가

"무시 가능한 수준(|ρ| < 0.2) 이라고 결론지을 수 있다"

── 반면 n = 50 이었다면 ──
r = 0.04,  95% CI [-0.24, 0.31]
등가 검정: p = 0.14  → 결론 유보
"상관이 없다고 말할 근거가 없다"

import numpy as np
from scipy import stats

za = stats.norm.ppf(0.975)
print("H0: ρ=0, 양측 0.05")
print(f"{'ρ':>6s} {'검정력 0.80':>11s} {'0.90':>8s} {'0.95':>8s} {'0.99':>8s}")
for rho in [0.1, 0.2, 0.3, 0.5, 0.7, 0.9]:
    row = []
    for pw in [0.80, 0.90, 0.95, 0.99]:
        zb = stats.norm.ppf(pw)
        row.append(int(np.ceil(((za + zb) / np.arctanh(rho))**2 + 3)))
    print(f"{rho:6.1f} {row[0]:11d} {row[1]:8d} {row[2]:8d} {row[3]:8d}")

print("\n주어진 n 에서 검출 가능한 최소 상관 (검정력 0.80)")
zb = stats.norm.ppf(0.80)
print(f"{'n':>6s} {'최소 ρ':>9s}")
for n in [20, 30, 50, 100, 200, 500, 1000]:
    print(f"{n:6d} {np.tanh((za + zb) / np.sqrt(n - 3)):9.4f}")

H0: ρ=0, 양측 0.05
     ρ    검정력 0.80     0.90     0.95     0.99
   0.1         783     1047     1294     1828
   0.2         194      259      320      451
   0.3          85      113      139      195
   0.5          30       38       47       64
   0.7          14       17       21       28
   0.9           7        8        9       12

주어진 n 에서 검출 가능한 최소 상관 (검정력 0.80)
     n      최소 ρ
    20    0.5912
    30    0.4924
    50    0.3873
   100    0.2770
   200    0.1970
   500    0.1250
  1000    0.0885

rng = np.random.default_rng(31001)
B = 20_000
print("유의한 연구만 보면 r 이 얼마나 부풀려지나")
print(f"{'참 ρ':>6s} {'n':>5s} {'유의 확률':>9s} "
      f"{'유의한 경우 평균 |r|':>19s} {'부풀림 배율':>11s}")
for rho, n in [(0.2, 40), (0.2, 100), (0.2, 200), (0.1, 40), (0.5, 20)]:
    rs, sig = [], 0
    for _ in range(B):
        z = rng.standard_normal((n, 2))
        r = np.corrcoef(z[:, 0], rho * z[:, 0]
                        + np.sqrt(1 - rho**2) * z[:, 1])[0, 1]
        t = r * np.sqrt((n - 2) / (1 - r**2))
        if 2 * stats.t.sf(abs(t), n - 2) < 0.05:
            sig += 1
            rs.append(abs(r))
    print(f"{rho:6.1f} {n:5d} {sig / B:9.4f} {np.mean(rs):19.4f} "
          f"{np.mean(rs) / rho:11.2f}")

유의한 연구만 보면 r 이 얼마나 부풀려지나
   참 ρ     n     유의 확률       유의한 경우 평균 |r|      부풀림 배율
   0.2    40    0.2351              0.3964        1.98
   0.2   100    0.5155              0.2736        1.37
   0.2   200    0.8139              0.2222        1.11
   0.1    40    0.0927              0.3783        3.78
   0.5    20    0.6382              0.5982        1.20

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(30004)
B = 5_000
print("두 변수가 독립인데 군집 효과가 있는 자료 (전체 n=200)")
print(f"{'군집 수':>7s} {'군집 크기':>9s} {'ICC':>6s} {'1종 오류':>9s} "
      f"{'설계효과':>9s} {'유효 n':>8s}")
for k, m in [(200, 1), (50, 4), (20, 10), (10, 20), (5, 40)]:
    for icc in ([0.0] if m == 1 else [0.3]):
        a = 0
        for _ in range(B):
            gx = rng.standard_normal(k) * np.sqrt(icc)
            gy = rng.standard_normal(k) * np.sqrt(icc)
            x = np.repeat(gx, m) + rng.standard_normal(k * m) * np.sqrt(1 - icc)
            y = np.repeat(gy, m) + rng.standard_normal(k * m) * np.sqrt(1 - icc)
            a += stats.pearsonr(x, y).pvalue < 0.05
        deff = 1 + (m - 1) * icc
        print(f"{k:7d} {m:9d} {icc:6.1f} {a / B:9.4f} {deff:9.2f} "
              f"{200 / deff:8.1f}")

두 변수가 독립인데 군집 효과가 있는 자료 (전체 n=200)
   군집 수     군집 크기    ICC     1종 오류      설계효과     유효 n
    200         1    0.0    0.0478      1.00    200.0
     50         4    0.3    0.0820      1.90    105.3
     20        10    0.3    0.1510      3.70     54.1
     10        20    0.3    0.2262      6.70     29.9
      5        40    0.3    0.3024     12.70     15.7

학급 20개, 학급당 10명

  학급 평균 키 와 학급 평균 성적이 우연히 상관되면
  → 학생 200명 수준에서 "유의한 상관" 이 나온다
  → 실제로는 학급 20개의 우연일 뿐

학급 · 학교 · 병원 · 의사 · 가구 · 마을 · 회사
반복측정(같은 사람의 여러 관측)
실험의 배치(batch) · 조사원 · 실험일

import warnings
warnings.filterwarnings("ignore")

import numpy as np
from scipy import stats

rng = np.random.default_rng(30005)
B = 20_000
print("1종 오류 (명목 0.05, 이변량 정규)")
print(f"{'참 ρ0':>7s} {'n':>5s} {'피셔 z':>9s} {'t 공식을 오용':>13s}")
for rho0 in [0.0, 0.5, 0.9]:
    for n in [20, 100]:
        a = b = 0
        for _ in range(B):
            z = rng.standard_normal((n, 2))
            r = np.corrcoef(z[:, 0], rho0 * z[:, 0]
                            + np.sqrt(1 - rho0**2) * z[:, 1])[0, 1]
            zz = (np.arctanh(r) - np.arctanh(rho0)) * np.sqrt(n - 3)
            a += 2 * stats.norm.sf(abs(zz)) < 0.05
            t = (r - rho0) * np.sqrt(n - 2) / np.sqrt(1 - r**2)
            b += 2 * stats.t.sf(abs(t), n - 2) < 0.05
        print(f"{rho0:7.1f} {n:5d} {a / B:9.4f} {b / B:13.4f}")

1종 오류 (명목 0.05, 이변량 정규)
   참 ρ0     n      피셔 z      t 공식을 오용
    0.0    20    0.0485        0.0477
    0.0   100    0.0515        0.0513
    0.5    20    0.0502        0.0260
    0.5   100    0.0526        0.0244
    0.9    20    0.0517        0.0001
    0.9   100    0.0498        0.0000

ρ = 0   : r 의 분포가 0 을 중심으로 대칭  →  t 근사가 정확
ρ = 0.9 : r 의 분포가 왼쪽으로 크게 치우침 →  대칭 근사가 실패
          (r 이 1 을 넘을 수 없으므로 위쪽이 눌린다)

H0: ρ = 0 의 검정       →  분포 가정에 매우 견고
구간 · ρ0 ≠ 0 의 검정   →  이변량 정규성이 필요

"비정규라서 검정을 못 한다"는 대개 과장이고
"비정규라서 구간을 믿을 수 없다"는 대개 맞다

□ 어느 상관계수인가 (피어슨/스피어만/켄들)
□ n (결측 처리 후)
□ 신뢰구간 (방법 명시: 피셔 z / 부트스트랩)
□ 관측의 독립성 (군집·시계열 여부)
□ 검정력 또는 검출 가능한 최소 ρ
□ 여러 검정을 했다면 보정
□ 산점도

운동 빈도와 우울 점수 (n = 312)

  피어슨 r = -0.28,  95% CI [-0.38, -0.17],  t(310) = -5.14, p < 0.001
  스피어만 r_s = -0.26  (두 값이 가까워 이상점 영향은 작다)

  구간은 피셔 z 변환으로 구했다. 이변량 정규성은 QQ 그림으로
  점검했으며 뚜렷한 이탈은 없었다.

  이 표본 크기는 |ρ| ≥ 0.16 을 검정력 0.80 으로 검출한다.

  참가자는 12개 지역센터에서 모집되었으나 센터 수준의
  군집내상관은 0.01 미만이었다 (설계효과 1.2).