기상·기후 예측의 평가

먼저 읽을 장 27.1 기상·기후 데이터의 기초21.2 예측분포와 구간의 평가

기상학은 예측을 매일 검증받아 온 가장 오래된 분야 중 하나다. 내일의 날씨는 내일 확인되고, 틀린 예보는 기록에 남는다. 그 과정에서 기상 분야는 "좋은 예보"를 재는 지표와 그 지표의 함정을 함께 쌓아 왔다. 데이터 기반 모델이 이 분야에 들어오면서 그 함정들이 다시 중요해졌다. 평균 제곱 오차로 학습한 모델은 평균 제곱 오차로 평가하면 유리하고, 그 유리함은 흐릿한 예보라는 형태로 나타난다.

이 장은 기상·기후 예측을 평가하는 표준 도구를 한 표에 모으는 법을 다룬다. 결정론적 지표(위도 가중 RMSE와 ACC, 27.2장의 모델 비교에서도 쓴다), 확률적 지표(21.2장의 CRPS와 보정 진단을 앙상블에 적용), 에너지 스펙트럼으로 흐릿함을 드러내는 법, 극한값과 추세처럼 표본이 적은 통계의 불확실성, 그리고 물리 법칙과의 일관성이다. 적정 점수 규칙과 보정의 일반론은 21.2장에서 보았으므로, 여기서는 그것이 격자장과 앙상블과 긴 기후 기록에서 어떻게 달라지는지에 집중한다.

학습 목표

  • 기상 분야의 표준 지표를 사용한다
  • 흐릿한 예측의 함정을 피한다
  • 추세와 극한의 불확실성을 자기상관과 비정상성을 고려해 평가한다

결정론적 지표

왜 필요한가. 결정론적 예보는 하나의 장을 낸다. 그 장이 정답과 얼마나 다른지를 재는 방법은 여러 가지이고, 방법마다 보는 것이 다르다. 크기를 보는 지표와 패턴을 보는 지표를 함께 써야 예보의 성격이 드러난다.

위도 가중 RMSE. 27.1장에서 본 것처럼 위경도 격자에서는 칸의 면적이 cos⁡φ\cos\varphi에 비례한다.

RMSE=1N∑n=1N∑i,jwi (fnij−onij)2∑i,jwi,wi=cos⁡φi\mathrm{RMSE} = \sqrt{\frac{1}{N}\sum_{n=1}^{N}\frac{\sum_{i,j} w_i\,(f_{nij} - o_{nij})^2}{\sum_{i,j} w_i}}, \qquad w_i = \cos\varphi_i

ff는 예보, oo는 정답(분석장), nn은 예보 사례, i,ji, j는 위도·경도 색인이다. 사례마다 MSE를 구해 평균한 뒤 제곱근을 취하는 것이 흔하지만, 사례마다 RMSE를 구해 평균하는 문헌도 있어 숫자가 조금 달라진다. 정의를 함께 적는다. 편향(bias, 평균 오차)도 함께 보고한다. RMSE가 같아도 체계적 편향인지 무작위 오차인지는 다르다.

Anomaly correlation coefficient. ACC는 예보와 정답의 기후값 편차 사이의 위도 가중 공간 상관이다(27.2장).

ACC=∑i,jwi fij′ oij′∑i,jwi fij′2  ∑i,jwi oij′2,f′=f−c, o′=o−c\mathrm{ACC} = \frac{\sum_{i,j} w_i\, f'_{ij}\, o'_{ij}}{\sqrt{\sum_{i,j} w_i\, f'^2_{ij}\;\sum_{i,j} w_i\, o'^2_{ij}}}, \qquad f' = f - c,\ o' = o - c

cc는 학습 기간의 기후값이다. 기후값을 빼지 않고 원래 값의 상관을 내면 계절 주기와 남북 기온 차이 같은 "누구나 맞히는" 구조 때문에 상관이 0.99에 가깝게 나온다. 편차의 상관이어야 예보의 정보를 잰다. 상관이므로 예보 편차에 양의 상수를 곱해도 ACC는 변하지 않는다. ACC는 패턴을, RMSE는 크기와 패턴을 함께 본다.

두 지표를 같은 정답에 대한 두 예보로 비교한다. 위도권 한 바퀴(128칸)의 장 500개를 만든다. 동서 파수 kk의 진폭은 k−1.5k^{-1.5}(에너지 스펙트럼 ∝k−3\propto k^{-3})이고, 파수마다 예측 가능한 정도 ρk=e−k/6\rho_k = e^{-k/6}가 다르다. 큰 규모는 잘 맞고 작은 규모는 거의 맞지 않는다. 날카로운 예보는 예측 가능한 부분에 예측할 수 없는 부분의 다른 실현을 더한, 실제로 있을 법한 장이다. 흐릿한 예보는 예측 가능한 부분만 남긴 조건부 평균, 즉 MSE를 최소화하는 예보다.

import numpy as np

rng = np.random.default_rng(0)
N, L = 500, 128                                   # 사례 수, 위도권 한 바퀴의 격자점 수
k = np.arange(1, L // 2)                          # (63,) 동서 파수
amp = k ** -1.5                                   # 에너지 스펙트럼 E(k) ∝ k^-3
rho = np.exp(-k / 6)                              # 파수별 예측 가능성: 큰 규모일수록 잘 맞는다

def field(coef):                                  # coef: (N, 63) 복소 계수 → (N, L) 실수 장
    full = np.zeros((coef.shape[0], L // 2 + 1), complex); full[:, 1:L // 2] = coef
    return np.fft.irfft(full, n=L) * L

cplx = lambda: amp * (rng.normal(size=(N, k.size)) + 1j * rng.normal(size=(N, k.size))) / np.sqrt(2)
c_a, c_b, c_b2 = cplx(), cplx(), cplx()           # 예측 가능한 부분, 예측할 수 없는 부분과 그 다른 실현
truth = field(rho * c_a + np.sqrt(1 - rho**2) * c_b)
sharp = field(rho * c_a + np.sqrt(1 - rho**2) * c_b2)    # 실현 가능한 장 하나(스펙트럼이 정답과 같다)
blurry = field(rho * c_a)                         # 조건부 평균(MSE 최적): 작은 규모가 사라진다

acc = lambda f: np.mean((f * truth).sum(1) / np.sqrt((f**2).sum(1) * (truth**2).sum(1)))
for name, f in [("날카로운 예보", sharp), ("흐릿한 예보", blurry)]:
    print(f"{name}: RMSE {np.sqrt(((f - truth)**2).mean()):.3f}  ACC {acc(f):.3f}  표준편차 {f.std():.3f}")
print(f"정답 표준편차 {truth.std():.3f}")
날카로운 예보: RMSE 1.283  ACC 0.548  표준편차 1.553
흐릿한 예보: RMSE 0.900  ACC 0.725  표준편차 1.221
정답 표준편차 1.502

흐릿한 예보가 RMSE와 ACC 모두에서 이긴다. 날카로운 예보의 RMSE는 흐릿한 예보의 약 2\sqrt{2}배다. 예측할 수 없는 작은 규모를 날카로운 예보는 엉뚱한 위치에 그리고(오차 한 번), 실제 위치에는 그리지 않는다(오차 또 한 번). 이를 이중 벌점(double penalty)이라 한다. 흐릿한 예보는 그 규모를 아예 그리지 않아 벌점을 한 번만 받는다. 대신 흐릿한 예보의 표준편차(1.221)는 정답(1.502)보다 작다. 이 차이가 어디서 오는지는 "흐릿함의 함정" 절에서 스펙트럼으로 본다.

연구에서의 의미. 결정론적 지표는 변수(500 hPa 지위고도, 850 hPa 기온, 2 m 기온, 강수 등)와 리드 타임별로 표를 만든다. 데이터 기반 모델과 수치 예보를 비교한 논문의 "몇 %의 목표에서 더 낫다"는 요약은 이 표의 칸을 센 것이다. 어떤 변수·리드·정답을 넣었는지에 따라 그 비율은 달라진다.

확률적 지표

왜 필요한가. 앙상블은 예측분포를 표본으로 낸다. 앙상블 평균의 RMSE만 보면 분포의 정보를 버리고, 멤버 하나의 RMSE를 보면 앙상블이 불리해진다(27.2장). 분포 전체를 평가하는 지표와 분포의 보정(calibration)을 진단하는 도구가 따로 필요하다.

CRPS. 21.2장의 표본 CRPS를 격자점마다 계산하고 위도 가중으로 평균한다.

CRPS=1M∑m∣xm−y∣−12M2∑m,m′∣xm−xm′∣\mathrm{CRPS} = \frac{1}{M}\sum_{m}\lvert x_m - y\rvert - \frac{1}{2M^2}\sum_{m,m'}\lvert x_m - x_{m'}\rvert

xmx_m은 멤버 mm의 값, yy는 정답이다. 멤버가 하나면 둘째 항이 0이 되어 절댓값 오차가 된다. 그래서 결정론적 예보와 앙상블을 같은 척도에서 비교할 수 있다. 멤버 수가 적으면 이 추정량은 무한 앙상블의 CRPS보다 조금 크게 나오므로, 멤버 수가 다른 앙상블끼리 비교할 때는 둘째 항의 분모를 M(M−1)M(M-1)로 바꾼 "fair" 판을 쓰기도 한다.

Spread–skill 관계. 보정된 앙상블에서는 정답이 멤버들과 같은 분포에서 나온 것처럼 행동한다. 그러면 앙상블 평균의 오차 분산은 멤버 분산의 (M+1)/M(M+1)/M배가 된다. 정답과 평균의 차이 y−xˉy - \bar{x}의 분산이 σ2+σ2/M\sigma^2 + \sigma^2/M이기 때문이다. 그래서

RMSE(xˉ)≈M+1M  s‾,s‾=멤버 분산(불편)의 평균\mathrm{RMSE}(\bar{x}) \approx \sqrt{\tfrac{M+1}{M}}\;\overline{s}, \qquad \overline{s} = \sqrt{\text{멤버 분산(불편)의 평균}}

이 성립해야 한다(Fortin 외, 2014). 유한 앙상블 보정을 곱한 spread와 앙상블 평균 RMSE의 비, 즉 (M+1)/M s‾ / RMSE(xˉ)\sqrt{(M+1)/M}\,\overline{s}\,/\,\mathrm{RMSE}(\bar{x})가 1보다 작으면 앙상블이 과신(퍼짐 부족), 크면 과소신이다.

Rank histogram. 각 사례·격자점에서 정답이 MM개 멤버 사이의 몇 번째 순위에 들어가는지 세면 M+1M+1개 칸의 히스토그램이 된다. 보정된 앙상블이면 평평하다. U자이면 정답이 멤버 범위 밖에 너무 자주 떨어지는 퍼짐 부족, 가운데가 솟으면 퍼짐 과다, 한쪽으로 기울면 편향이다. 다만 평평한 히스토그램이 보정을 보장하지는 않고, U자는 조건부 편향으로도 생길 수 있다(Hamill, 2001). 21.2장의 PIT 히스토그램을 앙상블에 옮긴 것이다.

M = 10
def ensemble(scale):                              # scale < 1이면 예측할 수 없는 부분을 과소 표현
    return np.stack([field(rho * c_a + scale * np.sqrt(1 - rho**2) * cplx()) for _ in range(M)])  # (M, N, L)

def crps(ens):                                    # 격자점별 표본 CRPS의 평균
    xs, w = np.sort(ens, 0), 2 * np.arange(1, len(ens) + 1) - len(ens) - 1
    pair = 2 * (w[:, None, None] * xs).sum(0) / len(ens) ** 2
    return (np.abs(ens - truth).mean(0) - pair / 2).mean()

print(f"흐릿한 결정론 예보 CRPS(=MAE) {np.abs(blurry - truth).mean():.3f}")
print(f"날카로운 단일 예보 CRPS(=MAE) {np.abs(sharp - truth).mean():.3f}")
for name, s in [("보정된 앙상블", 1.0), ("퍼짐 부족 앙상블", 0.5)]:
    ens = ensemble(s)
    rm = np.sqrt(((ens.mean(0) - truth) ** 2).mean())
    sp = np.sqrt(ens.var(0, ddof=1).mean() * (M + 1) / M)          # 유한 앙상블 보정
    rank = (ens < truth).sum(0).ravel()                             # 정답의 순위 0..M
    hist = np.bincount(rank, minlength=M + 1) / rank.size
    print(f"{name}: CRPS {crps(ens):.3f}  spread/RMSE {sp / rm:.2f}  "
          f"순위 양 끝 비율 {hist[0] + hist[-1]:.2f} (평평하면 {2 / (M + 1):.2f})")
흐릿한 결정론 예보 CRPS(=MAE) 0.716
날카로운 단일 예보 CRPS(=MAE) 1.027
보정된 앙상블: CRPS 0.558  spread/RMSE 1.00  순위 양 끝 비율 0.18 (평평하면 0.18)
퍼짐 부족 앙상블: CRPS 0.574  spread/RMSE 0.52  순위 양 끝 비율 0.46 (평평하면 0.18)

세 가지를 읽는다.

  • CRPS에서는 날카로운 앙상블이 흐릿한 결정론 예보를 이긴다(0.558 대 0.716). RMSE에서 졌던 날카로운 장들이 여러 개 모여 분포를 이루면, 예측할 수 없는 부분의 불확실성을 정직하게 표현한 것으로 보상받는다.
  • 보정된 앙상블은 spread/RMSE가 1.00, rank histogram의 양 끝 비율이 평평한 값(0.18)과 같다. 퍼짐을 절반으로 줄인 앙상블은 spread/RMSE 0.52, 양 끝 비율 0.46으로 과신이 분명히 드러난다.
  • 그런데 CRPS는 두 앙상블을 거의 구별하지 못한다(0.558 대 0.574). CRPS는 적정 점수라 보정이 나빠지면 기대값이 나빠지지만, 그 차이가 작을 수 있다. 보정은 spread–skill과 rank histogram으로 따로 진단한다.

흐릿함의 함정

왜 필요한가. 앞의 두 절에서 흐릿한 예보는 RMSE와 ACC를 모두 이겼다. MSE로 학습한 모델은 정확히 이 방향으로 간다. 리드 타임이 길어질수록 예측할 수 없는 규모가 커지므로, 그 규모를 지우는 것이 MSE를 줄이는 최선의 전략이 된다. 결정론적 지표만 보면 이 흐릿함은 결함이 아니라 성능으로 보인다.

정의와 진단. 흐릿함은 공간 규모별 에너지로 드러난다. 장을 동서 방향으로 Fourier 변환한 계수 f^k\hat f_k로 파수별 에너지 E(k)=∣f^k∣2‾E(k) = \overline{\lvert \hat f_k \rvert^2}를 정의한다(윗줄은 사례·위도 평균). 실현 가능한 대기 상태는 정답과 같은 모양의 스펙트럼을 가져야 한다. 예보와 정답의 에너지 비 Ef(k)/Eo(k)E_f(k)/E_o(k)를 파수별로 그리면, 1 아래로 떨어지는 파수에서 그 규모의 변동이 사라졌다는 뜻이다. 구면 전체에서는 구면 조화 함수로 같은 계산을 한다.

def power(f):                                     # (N, L) → (L//2+1,) 파수별 평균 에너지
    return (np.abs(np.fft.rfft(f, axis=1)) ** 2).mean(0) / L**2

ks = [1, 2, 4, 8, 16, 32]
print("파수 k          " + "".join(f"{kk:7d}" for kk in ks))
for name, f in [("날카로운 예보", sharp), ("흐릿한 예보", blurry), ("보정된 앙상블 평균", ensemble(1.0).mean(0))]:
    print(f"{name:10s}" + "".join(f"{v:7.2f}" for v in power(f)[ks] / power(truth)[ks]))
파수 k                1      2      4      8     16     32
날카로운 예보      1.10   0.91   0.95   0.93   0.88   1.09
흐릿한 예보       0.72   0.50   0.26   0.07   0.00   0.00
보정된 앙상블 평균   0.75   0.54   0.33   0.15   0.10   0.12

날카로운 예보의 에너지 비는 모든 파수에서 1 근처다(사례 수가 유한해서 생기는 0.9–1.1의 흔들림). 흐릿한 예보는 ρk2\rho_k^2를 따라 떨어져, 파수 8에서 7%, 16 이상에서 0이다. 그 규모의 변동을 전혀 그리지 않는다. 앙상블 평균도 흐릿하다. 10개 멤버의 평균이라 예측할 수 없는 부분의 에너지가 약 1/10로 줄어 작은 규모에서 0.1 근처에 머문다. 앙상블 평균은 결정론적 지표용이고, 실현 가능한 장은 멤버다.

연구에서의 의미. 데이터 기반 예보의 평가표에는 리드 타임별 에너지 스펙트럼(또는 정답 대비 에너지 비)을 넣는다. 스펙트럼이 리드에 따라 작은 규모부터 깎여 나가면, RMSE의 개선 일부는 흐릿함에서 온 것이다. 이것은 응용에서 바로 문제가 된다. 태풍 중심 기압, 국지 호우의 강도, 풍력 발전에 중요한 돌풍은 모두 작은 규모의 세기다. 흐릿한 예보는 이들을 체계적으로 과소평가한다. 다중 step 학습(27.2장)이나 손실 함수의 선택이 스펙트럼을 어떻게 바꾸는지 함께 보고하고, 생성모델 기반 예보는 CRPS와 스펙트럼을 함께 내는 것이 기본이다.

극한과 사건

왜 필요한가. 사회적 피해는 평균이 아니라 극한에서 온다. 극한은 정의상 드물어서 표본이 적고, 평균 지표에는 거의 기여하지 않는다. 평균 RMSE가 좋은 모델이 폭염과 극한 강수를 체계적으로 놓칠 수 있다.

사건 기반 지표. "일강수 50 mm 이상"처럼 문턱으로 사건을 정의하면 예·아니오의 분할표(contingency table)로 평가한다. 적중률(probability of detection), 오경보율(false alarm ratio), 임계성공지수(critical success index) 같은 지표가 여기서 나온다. 드문 사건에서는 "항상 아니오"라는 예보가 정확도(accuracy)로는 거의 완벽하므로, 정확도 대신 이런 지표와 확률 예보의 Brier score(21.2장)를 쓴다. 강수처럼 위치 오차가 큰 변수는 격자점 일치 대신 주변 영역 안에서 사건 비율을 비교하는 이웃 기반 지표(fractions skill score 등)로 이중 벌점을 피한다.

극값 분포와 재현 수준. 해마다의 최댓값(연 최댓값)은 한 해 안의 표본 수(블록 크기)가 충분히 크면 일반화 극값 분포(generalized extreme value, GEV)로 근사된다(Coles, 2001).

F(z)=exp⁡ ⁣{−[1+ξ z−μσ]−1/ξ}F(z) = \exp\!\left\{-\left[1 + \xi\,\frac{z - \mu}{\sigma}\right]^{-1/\xi}\right\}

μ\mu는 위치, σ>0\sigma > 0는 척도, ξ\xi는 모양 모수다. ξ>0\xi > 0이면 꼬리가 두꺼운 Fréchet형, ξ<0\xi < 0이면 상한이 있는 Weibull형, ξ→0\xi \to 0이면 Gumbel형이다. 해마다 넘을 확률이 1/T1/T인 값 zT=F−1(1−1/T)z_T = F^{-1}(1 - 1/T)를 TT년 재현 수준(return level)이라 한다. "100년 빈도 폭염"은 100년마다 한 번 온다는 뜻이 아니라 해마다 1%의 확률로 온다는 뜻이다.

여름 92일의 일 최고기온(자기상관 0.7의 편차)에서 연 최댓값을 만들고, 같은 기후에서 뽑은 30년 표본 세 개에 GEV를 맞춘다. 아주 긴 모의(20만 년)의 경험적 분위수를 "참" 재현 수준으로 쓴다. scipy의 genextreme은 모양 모수를 c=−ξc = -\xi로 정의하므로 부호를 바꿔 읽는다(집필 시점 예시. 설치한 버전의 문서에서 정의를 확인한다).

from scipy import stats

def annual_max(g, n_year, warming_c_per_yr=0.0):  # 여름 92일의 일 최고기온(AR(1) 편차)에서 해마다 최댓값
    a = np.zeros(n_year); mx = np.full(n_year, -np.inf)
    for _ in range(92):
        a = 0.7 * a + g.normal(0, 3.0 * np.sqrt(1 - 0.7**2), n_year)
        mx = np.maximum(mx, a)
    return 30.0 + mx + warming_c_per_yr * np.arange(n_year)    # (n_year,) [°C]

g = np.random.default_rng(1)
ref = annual_max(g, 200_000)                      # 아주 긴 정상 기후: '참' 재현 수준의 기준
print(f"참 100년 재현 수준 {np.quantile(ref, 1 - 1 / 100):.2f} °C")
for i in range(3):                                # 같은 기후의 서로 다른 30년 표본
    xi_neg, loc, scale = stats.genextreme.fit(annual_max(g, 30))   # scipy의 c는 -ξ
    rl = lambda T: stats.genextreme.ppf(1 - 1 / T, xi_neg, loc, scale)
    print(f"30년 표본 {i}: ξ = {-xi_neg:+.2f}  10년 {rl(10):.2f}  100년 {rl(100):.2f} °C")
level = np.quantile(ref, 1 - 1 / 100)
for warm in [0.5, 1.0, 2.0]:                      # 분포 전체가 평행이동한 미래 기후
    p = np.mean(ref + warm > level)
    print(f"평균 +{warm:.1f} °C: 과거 100년 수준을 넘을 확률 {p:.3f}/년 (재현 기간 {1 / p:.0f}년)")
참 100년 재현 수준 40.98 °C
30년 표본 0: ξ = -0.15  10년 38.74  100년 40.45 °C
30년 표본 1: ξ = -0.10  10년 38.50  100년 40.60 °C
30년 표본 2: ξ = +0.11  10년 39.71  100년 44.97 °C
평균 +0.5 °C: 과거 100년 수준을 넘을 확률 0.018/년 (재현 기간 54년)
평균 +1.0 °C: 과거 100년 수준을 넘을 확률 0.032/년 (재현 기간 31년)
평균 +2.0 °C: 과거 100년 수준을 넘을 확률 0.094/년 (재현 기간 11년)

30년으로 100년을 말하는 것은 외삽이다. 같은 기후에서 뽑은 30년 표본 세 개가 모양 모수 ξ\xi를 -0.15에서 +0.11까지 다르게 추정했다. 10년 재현 수준은 약 1.2 °C 폭 안에 모이지만, 100년 재현 수준은 40.45에서 44.97 °C까지 4.5 °C 넘게 갈린다. 꼬리 모양의 작은 차이가 자료 범위 밖에서 크게 벌어진다. 재현 수준에는 반드시 신뢰구간(부트스트랩, 프로파일 우도)을 붙인다.

분포 이동과 비정상성. 기후가 바뀌면 "과거의 100년 수준"의 의미가 바뀐다. 분포 전체가 1 °C만 평행이동해도 과거 100년 수준을 넘을 확률은 1%에서 3.2%로, 재현 기간은 31년으로 줄어든다. 2 °C면 11년이다. 평균의 작은 이동이 꼬리의 확률을 몇 배로 바꾼다. 그래서 추세가 있는 기록에 정상 GEV를 맞추면 미래 위험을 과소평가한다. 비정상 극값 모형은 위치 모수 μ\mu(필요하면 σ\sigma도)를 시간이나 전구 평균 기온의 함수로 두고 맞춘다(Coles, 2001). 21.4장의 분포 이동이 극값에서 나타나는 모습이다.

트렌드와 비정상성

왜 필요한가. 기후 평가의 많은 질문이 추세 질문이다. "이 지역의 여름 기온은 유의하게 올랐는가", "모델의 추세는 관측의 추세와 다른가". 기후 시계열은 자기상관이 있어서, 독립 표본을 가정한 표준오차는 추세를 실제보다 확실하게 보이게 한다.

정의. 시계열 yty_t에 최소제곱 직선을 맞춘 기울기 bb의 순진한 표준오차는 잔차가 독립이라고 가정한다. 잔차의 lag-1 자기상관을 r1r_1이라 하면, AR(1) 근사 아래 독립 표본으로 환산한 유효 표본 수는

neff=n 1−r11+r1n_\text{eff} = n\,\frac{1 - r_1}{1 + r_1}

이다. 표준오차의 자유도 n−2n - 2를 neff−2n_\text{eff} - 2로 바꿔 표준오차를 (n−2)/(neff−2)\sqrt{(n-2)/(n_\text{eff}-2)}배 키우고, tt 분포의 자유도도 neff−2n_\text{eff} - 2로 쓴다(Santer 외, 2000). 자기상관이 있는 두 시계열의 상관이나 공간장의 유효 자유도에도 같은 생각이 쓰인다(Bretherton 외, 1999). 11.7장의 AR(1) 과정과 16.7장의 자기상관 함수를 추세 검정에 적용한 것이다.

추세가 없는 AR(1) 연평균 편차(ϕ=0.6\phi = 0.6, 50년) 4000개에 두 검정을 적용해, 5% 유의수준에서 "유의한 추세"를 잘못 찾는 비율을 센다.

def trend_t(y):                                   # y: (R, n) → 기울기, 순진한 t, 보정한 t, n_eff
    n = y.shape[1]; t = np.arange(n) - (n - 1) / 2
    b = (y * t).sum(1) / (t**2).sum()
    res = y - y.mean(1, keepdims=True) - b[:, None] * t
    se = np.sqrt((res**2).sum(1) / (n - 2) / (t**2).sum())
    r1 = (res[:, 1:] * res[:, :-1]).sum(1) / (res**2).sum(1)          # 잔차의 lag-1 자기상관
    n_eff = np.clip(n * (1 - r1) / (1 + r1), 3, n)
    return b, b / se, b / (se * np.sqrt((n - 2) / (n_eff - 2))), n_eff

g = np.random.default_rng(4)
R, n, phi = 4000, 50, 0.6                         # 추세가 없는 AR(1) 연평균 편차 4000개, 각 50년
y = np.zeros((R, n)); y[:, 0] = g.normal(0, 1, R)
for t in range(1, n):
    y[:, t] = phi * y[:, t - 1] + np.sqrt(1 - phi**2) * g.normal(0, 1, R)
b, t_naive, t_adj, n_eff = trend_t(y)
crit = stats.t.ppf(0.975, n - 2)
print(f"유효 표본 수 중앙값 {np.median(n_eff):.1f} / {n}")
print(f"추세 없음인데 '유의'(5%): 순진한 검정 {np.mean(np.abs(t_naive) > crit):.3f}"
      f"  보정한 검정 {np.mean(np.abs(t_adj) > stats.t.ppf(0.975, n_eff - 2)):.3f}")
유효 표본 수 중앙값 16.4 / 50
추세 없음인데 '유의'(5%): 순진한 검정 0.320  보정한 검정 0.079

순진한 검정은 추세가 없는 자료의 32%에서 "유의한 추세"를 찾는다. 5%여야 할 거짓 양성이 6배다. 유효 표본 수로 보정하면 7.9%로 크게 줄지만 5%까지는 내려가지 않는다. 50년처럼 짧은 기록에서는 잔차로 추정한 r1r_1이 참값(0.6)보다 작게 나오는 경향이 있어 보정이 덜 되기 때문이다. 이론값 neff=50×0.4/1.6=12.5n_\text{eff} = 50 \times 0.4/1.6 = 12.5보다 추정한 중앙값 16.4가 큰 것이 그 흔적이다. 짧은 기록의 추세 검정은 보정한 뒤에도 보수적으로 읽는다.

비정상성과 평가 기간. 추세가 있으면 평가 지표도 기간에 따라 바뀐다. 기후값 baseline은 학습 기간의 평균이므로, 시험 기간이 더 따뜻하면 기후값 예보에 체계적 편향이 생기고, 그만큼 모델의 기후값 대비 스킬이 부풀려진다. 반대로 학습 기간보다 따뜻한 시험 기간에서 데이터 기반 모델이 처음 보는 극값을 만나면 성능이 떨어질 수 있다(21.4장). 평가 기간이 학습 기간과 기후적으로 얼마나 다른지를 함께 적는다.

연구에서의 의미. "모델이 관측된 추세를 재현한다"는 주장은 두 추세의 차이와 그 불확실성(자기상관 보정 포함)으로 검정한다. 기후 모델의 단일 모의는 내부 변동성 때문에 관측과 다른 추세를 낼 수 있으므로, 여러 초기 상태의 앙상블에서 추세의 분포를 보고 관측 추세가 그 안에 드는지 확인한다.

물리적 일관성

왜 필요한가. 데이터 기반 모델은 물리 법칙을 명시적으로 따르지 않는다. 지표가 좋은 예보도 질량이 새어 나가거나, 강수가 음수이거나, 바람과 기압의 관계가 깨져 있을 수 있다. 이런 오류는 짧은 리드의 RMSE에는 거의 보이지 않지만, 긴 rollout과 기후 규모의 적분에서 누적되어 결과를 망친다.

균형 관계. 중위도의 큰 규모 흐름에서는 바람이 등압선(또는 등지위고도선)을 따라 부는 지균 균형(geostrophic balance)이 근사적으로 성립하고, 연직 방향으로는 기압 경도력과 중력이 맞서는 정역학 균형(hydrostatic balance)이 성립한다. 예보한 지위고도장에서 계산한 지균풍과 예보한 바람을 비교하면, 두 변수가 서로 일관된 대기를 그리는지 볼 수 있다. 흐릿함은 이 관계도 바꾼다. 지위고도장이 평활되면 그 기울기로 계산한 바람이 약해진다.

보존량. 건조 공기의 전 지구 질량(지표 기압의 위도 가중 전구 평균에 비례), 전 지구 수분 수지(증발과 강수의 균형), 에너지는 짧은 기간에 거의 일정해야 한다. 27.1장의 위도 가중 평균으로 이 양들을 리드 타임별로 계산해, 정답의 변동 범위를 벗어나 표류(drift)하는지 확인한다. 수십 일 이상 굴리는 모델이라면 표류가 가장 먼저 보이는 진단이다.

물리적 범위. 강수·습도·해빙 농도는 음수가 될 수 없고, 상대습도와 해빙 농도는 상한이 있다. 이 범위를 벗어난 값의 비율을 표에 적고, 출력층에서 제약을 걸었다면 그 방식을 밝힌다. 24장의 물리 정보 모델은 이런 제약을 구조나 손실에 넣는 방법을 다룬다.

연구에서의 의미. 물리적 일관성 진단은 모델을 순위 매기는 지표라기보다 결과를 믿을 수 있는지를 거르는 검사다. 지표가 아무리 좋아도 보존량이 표류하는 모델의 장기 rollout으로 기후 통계를 내서는 안 된다.

실습

같은 예보에 RMSE와 에너지 스펙트럼·CRPS를 함께 적용해, RMSE만으로는 보이지 않는 차이를 찾는다.

  1. 예보 세 가지를 준비한다. 결정론적 모델 하나(27.2장의 실습에서 만든 모델이나 공개 벤치마크의 예보), 그 모델의 초기조건 앙상블(또는 생성모델 앙상블), 그리고 같은 정답에 대한 지속성·기후값 baseline. 공개 벤치마크의 예보 자료를 내려받아 써도 된다.
  2. 평가표를 만든다. 리드 타임별로 다음 열을 채운다. 위도 가중 RMSE·편향·ACC(결정론), CRPS·spread/RMSE·rank histogram의 양 끝 비율(확률), 몇 개 파수의 에너지 비(스펙트럼), 상위 1% 사건의 적중률·오경보율(극한), 전구 평균 지표 기압의 표류(물리).
  3. 흐릿함을 찾는다. RMSE 순위와 에너지 비 순위가 어긋나는 리드를 표시하고, 앙상블 평균과 단일 멤버의 스펙트럼을 비교한다.
  4. 추세와 극한. 관측(또는 재분석) 기록의 연평균 편차에서 추세를 구하고 자기상관을 보정한 신뢰구간을 붙인다. 연 최댓값에 GEV를 맞춰 재현 수준과 그 부트스트랩 구간을 낸다.

아래 두 연습문제가 2번과 4번의 핵심 계산이다.

실패 사례

  • RMSE만으로 흐릿한 예보를 더 좋다고 판단하는 것. 예측할 수 없는 규모를 지운 예보는 RMSE와 ACC를 모두 이긴다. 에너지 스펙트럼과 CRPS를 함께 보지 않으면 흐릿함이 성능으로 보고된다.
  • 앙상블을 CRPS 하나로만 평가하는 것. CRPS는 보정 오류에 둔감할 수 있다. spread–skill과 rank histogram을 함께 낸다.
  • 자기상관을 무시하고 추세의 유의성을 보고하는 것. 독립 표본을 가정한 검정은 추세가 없는 자료에서도 유의한 추세를 자주 찾는다.
  • 짧은 기록으로 긴 재현 수준을 점추정만 하는 것. 30년 자료의 100년 재현 수준은 외삽이고, 같은 기후의 다른 표본은 몇 °C 다른 값을 낸다.
  • 추세가 있는 기록에 정상 극값 분포를 맞추는 것. 과거 확률로 미래 위험을 말하게 된다.
  • 물리적 일관성을 확인하지 않고 긴 rollout을 기후 통계로 쓰는 것. 보존량의 표류는 짧은 리드의 지표에 보이지 않는다.

통과 기준

  • 결정론(위도 가중 RMSE·편향·ACC)·확률(CRPS·spread/RMSE·rank histogram)·스펙트럼(파수별 에너지 비)·극한(사건 지표 또는 재현 수준) 지표를 리드 타임별로 포함한 평가표를 만든다.
  • 평가표에서 RMSE 순위와 스펙트럼·CRPS 순위가 어긋나는 경우를 하나 이상 찾아 그 이유를 설명한다.
  • 연습문제 두 개(zonal-spectrum, trend-effective-n)의 테스트를 통과한다.

참고 문헌

  • Hamill (2001). Interpretation of rank histograms for verifying ensemble forecasts. Monthly Weather Review, 129(3), 550–560. — rank histogram의 해석과 함정
  • Fortin, Abaza, Anctil, Turcotte (2014). Why should ensemble spread match the RMSE of the ensemble mean? Journal of Hydrometeorology, 15, 1708–1713. — spread–skill 관계와 유한 앙상블 보정
  • Coles (2001). An Introduction to Statistical Modeling of Extreme Values. Springer. — GEV, 재현 수준, 비정상 극값 모형
  • Santer 외 (2000). Statistical significance of trends and trend differences in layer-average atmospheric temperature time series. Journal of Geophysical Research, 105(D6), 7337–7356. — 자기상관을 보정한 추세 검정
  • Bretherton 외 (1999). The effective number of spatial degrees of freedom of a time-varying field. Journal of Climate, 12(7), 1990–2009. — 유효 자유도