MCMC 추정과 적합도 검정

1. 오늘의 질문 (The Question)

5-6에서 우리는 아주 불편한 사실을 손으로 확인했다. triangle이나 gwesp가 들어오면 정규화 상수

가 쌍마다의 곱으로 쪼개지지 않는다. $n=4$짜리 망에서 64개를 전부 더해 $\kappa = 23.39495646$을 얻었을 때, 그 값이 6개 인수의 곱과 다르다는 것을 눈으로 봤다.

그런데 5-6에서 우리는 ergm(nwk ~ edges + gwesp(0.5, fixed=TRUE))를 실행해서 계수와 표준오차를 뽑아 썼다. 그러면 이상하다.

답은 순서대로 MCMC, mcmc.diagnostics(), gof()다. 오늘은 이 셋을 손으로 여섯 걸음 걸어 보는 것부터 시작한다.

2. 적률 방정식 다시 보기 (The Moment Equation Revisited)

5-1에서 유도했던 것을 다시 쓴다. 관측망 $y_{\text{obs}}$의 로그가능도는

이고, $\theta_k$로 미분할 때 두 번째 항은 이렇게 풀린다.

그래서 $\partial\ell/\partial\theta_k = 0$은 다음과 같다.

이것이 왜 중요한지는 §12에서 다시 나온다. 모형에 넣은 통계량은 정의상 맞을 수밖에 없다는 사실이 gof()를 읽는 핵심이기 때문이다.

3. 왜 그냥 못 푸는가 (Why It Cannot Be Solved Directly)

$E_{\boldsymbol\theta}[\mathbf g(Y)]$를 계산하려면 가능한 모든 망을 다 훑어야 한다. 무방향 망에서 쌍의 개수는 $N = \binom{n}{2}$이고, 각 쌍이 0 또는 1이니 망은 $2^N$개다.

그런데 자세히 보면 우리가 필요한 건 $\kappa$ 자체가 아니라 $E_{\boldsymbol\theta}[\mathbf g(Y)]$라는 평균이다. 평균은 다 더하지 않아도 구할 수 있다 — 표본을 뽑으면 된다. 전교생 1000명의 평균 키를 재려고 1000명을 다 재지 않고 무작위 50명을 재는 것과 같은 발상이다.

문제는 "$P_{\boldsymbol\theta}(y)$에서 망을 무작위로 뽑는 법"인데, 확률을 알려면 $\kappa$가 필요하고 $\kappa$를 모르니 뽑을 수가 없다 — 처럼 보인다. 이 고리를 끊는 것이 다음 절이다.

4. Metropolis–Hastings — $\kappa$가 약분된다 (Metropolis–Hastings: The Constant Cancels)

핵심 아이디어는 이렇다. 확률을 직접 계산하지 말고, 두 망의 확률 "비"만 쓰자. 비에서는 $\kappa$가 위아래로 똑같이 있어서 약분된다.

4.1 기호 (Notation)

5-3에서 쓴 표기를 그대로 가져온다. 쌍 $(i,j)$만 다르고 나머지는 같은 두 망을

4.2 비를 계산하면 $\kappa$가 사라진다 (The Ratio Kills the Constant)

4.3 걸음 규칙 (The Walk)

현재 망을 $y$라 하자. 한 걸음은 이렇게 걷는다.

③의 비를 §4.2에 넣으면 이렇게 갈린다.

4.4 왜 이게 맞는 답을 주는가 (Why This Works)

이 걸음은 세부 균형(detailed balance)을 만족한다. 즉 임의의 두 이웃 망 $y, y'$에 대해

제안 $q$는 대칭이다 — $y$에서 $(i,j)$를 골라 $y'$로 갈 확률도 $1/N$, 거꾸로도 $1/N$. 그래서 $q$가 양변에서 지워지고, 확인할 것은 $P(y)a(y\!\to\!y') = P(y')a(y'\!\to\!y)$뿐이다. $r = P(y')/P(y)$라 두고 $r \le 1$인 경우를 보면

양변이 같다. $r > 1$이면 역할만 바뀌므로 마찬가지다. 세부 균형이 성립하면 $P_{\boldsymbol\theta}$가 이 걸음의 정상분포(stationary distribution)이고, 충분히 오래 걸으면 방문 빈도가 $P_{\boldsymbol\theta}$에 수렴한다.

5. 손 계산 ① — 4명 망을 여섯 걸음 (Six Steps by Hand)

5-6 §5에서 전수 계산했던 바로 그 모형을 쓴다. $n=4$, 항은 edges + triangle,

변화 통계량은 5-6에서 구한 대로다.

즉 edges의 변화 통계량은 언제나 1이고, triangle의 변화 통계량은 $i$와 $j$의 공통 친구 수다.

따라서 켜는 제안의 지수는

여기서 $p_{ij}$는 그 간선을 뺀 상태에서 센 공통 친구 수다.

5.1 미리 계산해 둘 값 (Values to Prepare)

끄는 제안은 부호만 뒤집는다.

5.2 여섯 걸음 전개 (The Six Steps Expanded)

빈 망에서 출발한다. 쌍은 combn(4,2) 순서로 1=(1,2), 2=(1,3), 3=(1,4), 4=(2,3), 5=(2,4), 6=(3,4)이고, 각 걸음의 난수는 미리 뽑아 고정했다.

5.3 각 걸음을 말로 (Step by Step in Words)

6. R 검증 ① — 전수 계산을 재현하는가 (Verification I: Does the Walk Reproduce the Truth?)

$n=4$는 전수 계산이 가능하니, MCMC가 맞는 답으로 가는지 직접 채점할 수 있다. 5-6에서 64개를 다 더해 얻은 정확값은

## §5의 규칙을 그대로 코드로 옮긴 것
run_mh <- function(nstep, seed) {
  set.seed(seed); A <- matrix(0,4,4); Ee <- integer(nstep); Tt <- integer(nstep)
  for (s in 1:nstep) {
    k <- sample.int(6,1); i <- pr4[k,1]; j <- pr4[k,2]      # ① 쌍 하나 균등 선택
    cur <- A[i,j]; Am <- A; Am[i,j] <- Am[j,i] <- 0
    d_t <- sum(Am[i,]*Am[j,])                                # 공통 친구 수 p_ij
    logr <- if (cur==0) (th_e + th_t*d_t) else -(th_e + th_t*d_t)   # ③ 부호가 갈린다
    if (log(runif(1)) < logr) { A[i,j] <- A[j,i] <- 1-cur }  # ④ u < a 이면 뒤집기
    Ee[s] <- sum(A)/2; Tt[s] <- round(sum(diag(A%*%A%*%A))/6)
  }
  list(e=Ee, t=Tt)
}
r <- run_mh(100000, 20260803)

6.1 평균이 수렴하는가 (Convergence of the Mean)

6.2 분포 전체가 맞는가 (Does the Whole Distribution Match?)

평균만 맞추는 건 쉽다. 진짜 확인은 분포 전체다.

최대 오차가 0.00123이다. 확률 0.17짜리 칸에서 0.001 틀린 것이니 상대오차 0.7%. 전수 계산을 한 번도 하지 않고, 오직 "쌍 하나 뒤집고 지수 하나 계산하기"만 10만 번 반복해서 $2^6$개 망의 분포를 복원했다.

6.3 걸음은 서로 독립이 아니다 (The Steps Are Not Independent)

여기에 함정이 하나 있다. 한 걸음은 쌍 하나만 바꾸므로 연속한 표본이 거의 같은 망이다. 자기상관을 재 보면

바로 옆 걸음과는 0.83이나 붙어 있다. 6쌍짜리 망이라 6걸음쯤 지나야 겨우 풀린다. 그래서 10만 걸음이 10만 개의 독립 표본이 아니다.

7. MCMLE — 중요도 표집 (MCMLE: Importance Sampling)

이제 표본은 뽑을 수 있다. 그런데 우리가 원하는 건 $\boldsymbol\theta$를 찾는 것이고, $\ell(\boldsymbol\theta)$를 여러 $\boldsymbol\theta$에서 비교해야 한다. $\boldsymbol\theta$를 바꿀 때마다 MCMC를 새로 돌리면 너무 느리다.

Geyer–Thompson의 요령은 한 번 뽑은 표본을 재활용하는 것이다.

7.1 유도 (Derivation)

어떤 시작값 $\boldsymbol\theta_0$에서 표본을 뽑았다고 하자. 그러면

가운데 등호에서 한 일은 $e^{\boldsymbol\theta^\top g}$를 $e^{\boldsymbol\theta_0^\top g}\times e^{(\boldsymbol\theta-\boldsymbol\theta_0)^\top g}$로 쪼갠 것뿐이다. 그러면 앞부분이 정확히 $P_{\boldsymbol\theta_0}(y)$의 모양이 되고, 합이 기댓값으로 바뀐다. 이 기댓값은 $\boldsymbol\theta_0$에서 뽑은 표본으로 근사할 수 있다.

그리고 로그가능도의 차에서는 우리가 모르는 $\kappa(\boldsymbol\theta_0)$마저 지워진다.

7.2 왜 반복해야 하는가 (Why It Must Iterate)

8. R 검증 ② — 답을 아는 문제로 채점 (Verification II: Grading Against a Known Answer)

MCMC를 믿으려면 답을 아는 문제로 채점해 봐야 한다. 가라테에 edges만 넣은 모형이 그런 문제다. 5-2에서 손으로 풀었듯이 답이 닫힌 식으로 나온다.

8.1 중요도 표집이 $\kappa$ 비를 맞추는가 (Does Importance Sampling Get the Ratio Right?)

$\boldsymbol\theta_0 = -1.5$에서 2000개를 모의하고, §7.1 공식으로 $\log[\kappa(\theta)/\kappa(\theta_0)]$를 여러 $\theta$에서 추정해 정확값과 비교했다. edges만 있는 모형은 $\kappa$가 닫힌 식 $(1+e^{\theta})^{N}$이므로 정확값을 안다.

th0 <- -1.5
sim0 <- simulate(nwk ~ edges, coef=th0, nsim=2000, seed=7,
                 control=control.simulate(MCMC.burnin=10000, MCMC.interval=200))
g0 <- sapply(sim0, network.edgecount)
logIS <- function(th) { d <- (th-th0)*g0; mx <- max(d); mx + log(mean(exp(d-mx))) }
lhat  <- function(th) (th-th0)*mobs - logIS(th)

이 추정 곡선을 최대로 만드는 값은 $\hat\theta = -1.814238984$로, 정확한 MLE $-1.823307827$과 0.00906884 차이다. 표준오차 0.122의 7% 수준이니 실용적으로 무시할 만하다.

8.2 ergm에게 억지로 MCMC를 시켜 보면 (Forcing ergm to Use MCMC)

ergm은 edges만 있으면 똑똑하게 닫힌 식으로 풀어 버린다. force.main=TRUE로 MCMC 경로를 강제하면 채점이 된다.

ergm(nwk ~ edges, control=control.ergm(force.main=TRUE, MCMLE.maxit=20, seed=11))

9. mcmc.diagnostics() 읽는 법 (Reading the Diagnostics)

5-6에서 적합했던 edges + gwesp(0.5, fixed=TRUE)로 돌아가자. 이번엔 적합 과정을 들여다본다.

fit <- ergm(nwk ~ edges + gwesp(0.5, fixed=TRUE), control=control.ergm(seed=1))
mcmc.diagnostics(fit)

9.1 반복이 수렴하는 모습 (The Iteration Trace)

9.2 표본 궤적과 자기상관 (Trace and Autocorrelation)

mcmc.diagnostics()가 보고한 숫자는 이렇다.

10. 죽은 사슬은 어떻게 생겼나 (What a Dead Chain Looks Like)

5-6 §8에서 edges + triangle이 이렇게 죽었다.

Iteration 1 of at most 10:
  1 Optimizing with step length 0.4975.
  The log-likelihood improved by 6.6908.
  Estimating equations are not within tolerance region.
Iteration 2 of at most 10:
  Post-burnin sample is constant; returning.
Error: Unconstrained MCMC sampling did not mix at all. Optimization cannot continue.

이제 이 문장을 정확히 읽을 수 있다. "Post-burnin sample is constant"는 번인을 지난 뒤의 표본이 전부 같은 망이라는 뜻이다. 사슬이 한 망에 도착한 뒤 한 걸음도 못 움직였다.

어떤 망인지 직접 잡아 봤다. 5-6에서 구한 MPLE 계수 $\theta_{\text{edges}}=-2.6352326078$, $\theta_{\text{triangle}}=0.6877$로 관측망에서 출발해 2000걸음마다 하나씩 200개를 기록했다.

왜 못 나오는가? 완전망에서 간선 하나를 끄자는 제안을 생각해 보자. $n=34$에서 두 사람의 공통 친구는 나머지 32명 전부다.

11. gof()는 무엇을 채점하는가 (What GOF Grades)

계수가 나왔고 사슬도 건강하다고 하자. 그래도 남는 질문이 있다.

gof()의 절차는 간단하다.

④의 위치를 수치로 요약한 것이 MC $p$-value다. 관측값보다 크거나 같은 모의값의 비율과 작거나 같은 비율 중 작은 쪽을 두 배 한 값 (대략적으로는 "관측값이 얼마나 가장자리인가")이다.

11.1 무엇을 채점 대상으로 삼는가 (Which Features to Grade)

12. 모형 통계량이 맞는 것은 당연하다 (Model Statistics Are Matched by Construction)

§2의 적률 방정식을 다시 보자.

가라테 두 모형에서 실제로 확인해 보자.

13. R 검증 ③ — 가라테 두 모형 채점 (Verification III: Grading Two Karate Models)

gof(fit, GOF = ~ degree + espartners + distance + model,
    control = control.gof.ergm(nsim=200, seed=3))

13.1 공유 파트너 분포 — gwesp가 고친다 (ESP: Fixed by GWESP)

13.2 거리 분포 (Geodesic Distance Distribution)

edges만 있는 모형에서 거리 3인 쌍이 관측 137개인데 모의는 최소가 147이다. 200개 중 하나도 137 이하로 못 내려갔다. 뭉친 우정은 거리를 단축시키는데, 무작위 망은 그 단축을 만들지 못한 것이다.

13.3 차수 분포 — 둘 다 실패한다 (Degree: Both Models Fail)

14. 교실 적용 (Classroom Application)

14.1 발표 전에 반드시 확인할 것 (A Checklist Before Presenting)

14.2 각 gof 항을 교실 언어로 (GOF in Classroom Language)

15. 연습문제 (Exercises)

16. 해설과 답 (Solutions)

연습 1 해설 (Solution 1)

① 무엇을 보는가. 시작 망의 인접행렬은 이렇다. 색칠한 칸이 1인 자리다.

차수는 $(2,1,1,0)$, 간선 2개, 삼각형 0개다. 공통 친구 수 $p_{ij}$는 $i$행과 $j$행을 나란히 놓고 둘 다 1인 열을 세는 것이다 (단, 쌍 $(i,j)$ 자신은 0으로 지워 놓고 센다).

② 전개 — 0인 항도 모두 쓴다.

걸음 1: 쌍 (2,3), 현재 $A[2,3]=0$ → 켜는 제안

$u = 0.65 < 1$ → 수락. 망 = 1–2, 1–3, 2–3, 간선 3개, 삼각형 1개({1,2,3}).

걸음 2: 쌍 (1,2), 현재 $A[1,2]=1$ → 끄는 제안

먼저 $A[1,2]=A[2,1]=0$으로 지운 상태에서 공통 친구를 센다.

끄는 제안이므로 지수의 부호가 뒤집힌다.

$u = 0.91 > 0.7408182$ → 기각. 망은 그대로. 간선 3개, 삼각형 1개.

걸음 3: 쌍 (3,4), 현재 $A[3,4]=0$ → 켜는 제안

$u = 0.42 < 0.6065307$ → 수락. 망 = 1–2, 1–3, 2–3, 3–4, 간선 4개, 삼각형 1개.

연습 2 해설 (Solution 2)

① 무엇을 곱하는가. 계수 차부터 구한다.

삼각형 계수는 바뀌지 않았으므로 삼각형 수는 계산에 전혀 기여하지 않는다. 지수는 $-0.2 \times (\text{간선 수}) + 0 \times (\text{삼각형 수})$다.

② (가) 전개 — 0인 항도 모두.

③ (나)

④ (다) 관측망은 $\mathbf g(y_{\text{obs}}) = (3,\,1)$이므로

⑤ (라) 정확값과 비교 — 여기가 진짜 문제다.

연습 3 해설 (Solution 3)

① 표를 두 덩어리로 갈라 본다. 이것이 gof 표를 읽는 첫 동작이다.