Ch.11 Exercises — 반올림 오차의 누율구조·균등분포 누율·\(\delta\) 프로파일 이탈도 (McCullagh §11.7)

1 개요

Ch.11 의 본문(§11.2-§11.5)이 “모수가 선형 예측자 바깥에 있을 때” 의 추정을 다뤘다면, §11.7 의 네 연습문제는 그 주제에 직접 매달리지 않고 조금 비껴간 세 가지 부수 결과 를 묶어 놓는다. 겉보기엔 서로 무관해 보이지만 공통의 교훈은 하나다.

“유한한 정밀도로 기록된 관측치 는 이론적 확률 모형의 관측치와 다르다. 그 차이가 추론에 미치는 영향을 어디까지 계량할 수 있는가?”

문제	대응 본문 결과	핵심 메시지
11.1	Kolassa-McCullagh (1987) 반올림 정리	반올림 오차의 모든 누율을 간결하게 쓸 수 있다
11.2	감마 분산 함수의 반올림 보정	기록 정밀도가 분산 함수 선택에 영향
11.3	균등 분포의 네 누율	11.1 의 구체화 — 실제 계산에 필요한 상수
11.4	§11.4 선형화의 실전 적합	Ch.11 기법을 Ch.8 비선형 문제에 적용

처음 세 문제는 Ch.8 와 Ch.11 을 연결하는 이론적 접착제 다. §11.2 가 “분산 함수 안의 미지 모수” 를 다뤘다면, 11.2 연습문제는 “우리가 기록한 데이터의 분산 함수 자체가 이미 반올림 때문에 이상형(ideal form)에서 벗어나 있다” 는 사실을 보여 준다. 네 번째 문제는 §11.4 의 이중 IRLS 를 실제 데이터에 돌려서 프로파일 이탈도의 비이차성을 직접 확인시킨다.

2 문제 11.1 — 반올림 오차의 누율 구조

2.1 설정 — 무엇이 ‘작동’ 하는가

연속 확률변수 \(Z\) 의 밀도 \(f(z)\) 가 충분히 매끄럽다고 하자. 구체적으로 어떤 짝수 \(\nu \geq 2\) 에 대해

\[\int_{-\infty}^{\infty} |f^{(\nu)}(z)|\, dz < \infty\]

를 만족한다. 관측값이 정밀도 \(\epsilon = 10^{-d}\) 로 반올림되어 기록된다.

\[Y = \epsilon\,\lfloor Z/\epsilon \rceil, \qquad Z = Y + \epsilon R, \qquad R \in \left(-\tfrac{1}{2}, \tfrac{1}{2}\right]\]

여기서 \(\lfloor x \rceil\) 는 \(x\) 에 가장 가까운 정수, \(R\) 은 정규화 반올림 오차 다.

직관: \(Y, R\) vs \(Z\)

• \(Z\) 는 “이상적인 연속 관측치”. 실험실 용어로는 계측기가 무한 정밀도라면 얻을 값이다.
• \(Y\) 는 “실제 데이터” — 소수점 \(d\) 자리까지만 기록된 반올림 값.
• \(R = (Z-Y)/\epsilon\) 은 “우리가 버린 부분” 이다. \(R\) 은 \((-1/2, 1/2]\) 범위에 묶여 있다.

중요한 개념적 사실: \(R\) 은 \(Z\) 의 결정적 함수 다 (\(R = (Z/\epsilon) - \lfloor Z/\epsilon \rceil\)). 그러나 이 결정성에도 불구하고 \(\epsilon\) 이 작아질수록 \(R\) 과 \(Z\) 는 점근적으로 독립 이 된다. 이것이 Kolassa-McCullagh (1987) 의 놀라운 결과다.

2.2 결과 — 세 가지 주장

McCullagh-Nelder 는 Euler-Maclaurin 합산 공식으로 다음을 보이라고 요구한다.

1. \((Z, R)\) 의 결합 누율 모함수(CGF) 가 합 으로 분해된다: \[K_{Z,R}(\xi_1, \xi_2) = K_Z(\xi_1) + K_R(\xi_2) + O(\epsilon^\nu).\]

2. \(R\) 의 주변 CGF 가 주기적 함수의 로그로 주어진다: \[K_R(\xi) = \log\left(\frac{\sinh(\xi/2)}{\xi/2}\right) + O(\epsilon^\nu).\]

3. 반올림된 \(Y\) 의 누율은 홀수/짝수에 따라 다르게 변한다: \[\kappa_r(Y) \simeq \kappa_r(Z) \quad (r \text{ odd}), \qquad \kappa_r(Y) \simeq \kappa_r(Z) + \epsilon^r \kappa_r(R) \quad (r \text{ even}).\]

2.3 직관 — 왜 \(\sinh(\xi/2)/(\xi/2)\) 인가

\(R \sim U(-1/2, 1/2]\) 의 적률 모함수는 바로 쓸 수 있다.

\[ M_R(\xi) = E[e^{\xi R}] = \int_{-1/2}^{1/2} e^{\xi r}\, dr = \frac{e^{\xi/2} - e^{-\xi/2}}{\xi} = \frac{2\sinh(\xi/2)}{\xi}. \]

좌변에서 \(\xi = 0\) 일 때 1 이 되도록 다시 쓰면 \(M_R(\xi) = \sinh(\xi/2)/(\xi/2)\) — 여기에 로그를 취한 것이 바로 \(K_R(\xi)\) 다.

핵심은 “\(R\) 이 일반 \(U(-1/2, 1/2]\) 처럼 행동 한다” 는 것이다. 이상하다 — \(R\) 은 \(Z\) 의 결정적 함수인데 어떻게 \(Z\) 와 독립인 새 균등변수처럼 취급될 수 있는가? 답은 Euler-Maclaurin 의 핵심 메커니즘에 있다.

직관: Euler-Maclaurin 이 하는 일

합 \(\sum_{k} g(\epsilon k)\) 은 \(g\) 가 충분히 매끄러우면 적분 \((1/\epsilon)\int g(z)\, dz\) 와 지수적으로 작은 차이 밖에 나지 않는다 (구체적으로는 \(O(\epsilon^\nu)\)). 반올림을 CGF 의 푸리에 계수 언어로 옮기면, \(\epsilon\) 격자 위의 주기적 나머지 항 만 남는데, 이 항이 바로 \(K_R\) 의 “지역적” 기여분이다. 밀도 \(f\) 가 매끄러울수록 주기 항은 더 약해지고, 남는 것은 자율적으로 균등분포처럼 움직이는 \(R\) 의 기여뿐이다.

2.4 왜 홀수 누율은 변하지 않는가

\(R \sim U(-1/2, 1/2]\) 는 0 중심 대칭 이다. 대칭분포의 모든 홀수 누율은 0 이다:

\[\kappa_r(R) = 0 \quad \text{for } r \text{ odd}.\]

따라서 주장 3 에서 홀수 항은 \(\epsilon^r \cdot 0 = 0\) — \(Y\) 의 홀수 누율은 \(Z\) 의 그것과 같다. 평균(\(\kappa_1\)), 왜도 분자(\(\kappa_3\)) 등은 반올림해도 보존 된다.

짝수 누율은 다르다. 분산, 첨도 같은 짝수 누율은 \(\epsilon^r \kappa_r(R)\) 만큼 더해진다. \(\kappa_r(R)\) 은 11.3 에서 구체적으로 계산하게 되는 상수들이다.

2.5 \(\text{cov}(R, Y)\) 가 여전히 작지만 0 이 아닌 이유

주장 1 은 \(Z\) 와 \(R\) 이 독립 이라 말한다. 그러나 \(Y = Z - \epsilon R\) 이므로

\[\text{cov}(R, Y) = \text{cov}(R, Z) - \epsilon\,\text{var}(R) \approx 0 - \epsilon\,\text{var}(R) = -\epsilon/12.\]

(부호는 교재의 \(Y \approx Z - \epsilon R\) 표현에 맞춰 음수지만 원문은 양수 \(\epsilon/12\) 로 적혀있다 — 정의 \(Z = Y + \epsilon R\) 을 사용하면 \(\text{cov}(R, Y) = -\epsilon/12\).)

핵심 메시지는 부호가 아니라 크기 다. \(R\) 과 \(Z\) 는 \(O(\epsilon^\nu)\) 까지 독립이지만, \(R\) 과 \(Y\) 의 상관은 \(O(\epsilon)\) 수준으로 느슨하다. 수학적으로는 뉘앙스 차이 같지만, 실무에서는 “\(Y\) 의 분산과 상관 구조에서 반올림 성분 \(\epsilon R\) 이 1차 항으로 살아남는다” 는 뜻이다.

3 문제 11.2 — 감마 분포의 반올림 분산

3.1 유도

\(Z \sim \text{Gamma}(\nu, \mu/\nu)\) 로 평균 \(\mu_Z\), 분산 \(\mu_Z^2/\nu\) 이다. 관측치를 \(d\) 자릿수로 반올림한 \(Y\) 의 분산은 문제 11.1 의 결과로

\[ \text{var}(Y) = \kappa_2(Y) \simeq \kappa_2(Z) + \epsilon^2 \kappa_2(R) = \frac{\mu_Z^2}{\nu} + \frac{\epsilon^2}{12}. \]

(\(\kappa_2(R) = \text{var}(U) = 1/12\) 는 문제 11.3 에서 유도한다.)

3.2 왜 이것이 중요한가 — 분산 함수의 ‘바닥’

GLM 에서 감마 분산 함수는 \(V(\mu) = \mu^2\) 이다 — 평균이 작아지면 분산도 작아진다. 그런데 반올림 보정 항 \(\epsilon^2/12\) 는 \(\mu\) 에 무관한 상수 다. 결과적으로 \(Y\) 의 실제 분산 함수는

\[V_Y(\mu) = \frac{\mu^2}{\nu} + \frac{\epsilon^2}{12},\]

즉 \(\mu\) 가 \(\epsilon \sqrt{12/\nu}\) 보다 작은 영역에서는 바닥(floor) 이 생긴다. 이 영역에서는 감마 모형이 더 이상 맞지 않는다.

실무 함의: 정밀도가 분산 함수를 바꾼다

• 수명 1.5 시간 ±0.1 시간 정밀도로 기록된 소자 수명 데이터: \(\epsilon^2/12 = 8.3 \times 10^{-4}\). \(\nu = 10\) 이고 \(\mu = 1.5\) 일 때 \(\mu^2/\nu = 0.225\) — 반올림 기여는 \(0.4\%\) 로 무시 가능.
• 짧은 반응 시간 \(\mu = 0.01\) 초 를 \(10^{-2}\) 초 단위로 기록: \(\mu^2/\nu = 10^{-5}/\nu\) 인데 반올림 바닥은 \(8.3 \times 10^{-6}\) — 정밀도가 분산 함수를 지배 한다. 이때 감마 GLM 은 부적절하며 정규 근사가 더 적합할 수도 있다.

교재는 “반올림 바닥” 을 명시적으로 설계에 반영할 것을 권한다. 실측 데이터의 기록 정밀도 \(\epsilon\) 를 적어두고 분산 함수 검정 시 이 상수를 기준선으로 삼는다.

3.3 유효 자릿수로 반올림할 때

문제의 후속부: “\(Z\) 를 \(d\) 자릿수 유효 숫자(significant digits) 로 반올림하면 분산 함수가 어떻게 바뀌는가?”

유효 숫자 반올림은 \(\epsilon\) 가 \(\mu\) 에 비례 한다 — \(\epsilon \approx \mu \cdot 10^{-d}\). 따라서

\[ \text{var}(Y) \simeq \frac{\mu^2}{\nu} + \frac{(\mu \cdot 10^{-d})^2}{12} = \mu^2 \left(\frac{1}{\nu} + \frac{10^{-2d}}{12}\right). \]

결과는 여전히 \(V(\mu) = c \cdot \mu^2\) 형태 — 즉 감마 분산 함수를 보존한다. 단, 유효 “형태 모수” 가 \(\nu\) 에서

\[\nu_{\text{eff}} = \frac{\nu}{1 + \nu \cdot 10^{-2d}/12}\]

로 감소한다. \(d\) 자리 유효 숫자는 \(\nu\) 를 감소 시키지만 “바닥” 을 만들지는 않는다. 정규성 근사 관점에서는 유효 숫자 반올림이 더 안전 하다.

반올림 유형	분산 함수 변화	GLM 영향
\(d\) 소수점	\(\mu^2/\nu + \epsilon^2/12\) (상수 바닥)	작은 \(\mu\) 에서 분산 함수 형태 변경
\(d\) 유효 숫자	\(\mu^2(1/\nu + 10^{-2d}/12)\)	형태 보존, \(\nu\) 만 축소

4 문제 11.3 — 균등 분포의 네 누율

4.1 계산

\(U \sim U(0, 1)\) 의 적률은 직접 적분으로 쉽게 나온다.

\[\mu_1 = E[U] = \int_0^1 u\, du = \tfrac{1}{2},\] \[\mu_2 = E[U^2] = \tfrac{1}{3}, \quad \mu_3 = \tfrac{1}{4}, \quad \mu_4 = \tfrac{1}{5}.\]

누율은 적률의 관계식

\[\kappa_1 = \mu_1, \qquad \kappa_2 = \mu_2 - \mu_1^2,\] \[\kappa_3 = \mu_3 - 3\mu_1 \mu_2 + 2\mu_1^3,\] \[\kappa_4 = \mu_4 - 4\mu_1 \mu_3 - 3\mu_2^2 + 12\mu_1^2 \mu_2 - 6\mu_1^4\]

에서 계산된다.

\[ \kappa_1 = \tfrac{1}{2}, \qquad \kappa_2 = \tfrac{1}{3} - \tfrac{1}{4} = \tfrac{1}{12}, \]

\[ \kappa_3 = \tfrac{1}{4} - 3 \cdot \tfrac{1}{2} \cdot \tfrac{1}{3} + 2 \cdot \tfrac{1}{8} = \tfrac{1}{4} - \tfrac{1}{2} + \tfrac{1}{4} = 0, \]

\[ \kappa_4 = \tfrac{1}{5} - 4 \cdot \tfrac{1}{2} \cdot \tfrac{1}{4} - 3 \cdot \tfrac{1}{9} + 12 \cdot \tfrac{1}{4} \cdot \tfrac{1}{3} - 6 \cdot \tfrac{1}{16} = \tfrac{1}{5} - \tfrac{1}{2} - \tfrac{1}{3} + 1 - \tfrac{3}{8} = -\tfrac{1}{120}. \]

정리하면 \(U(0,1)\) 의 처음 네 누율은

\[\boxed{\kappa_1 = \tfrac{1}{2}, \quad \kappa_2 = \tfrac{1}{12}, \quad \kappa_3 = 0, \quad \kappa_4 = -\tfrac{1}{120}.}\]

4.2 CGF 로부터 교차 검증

앞서 \(M_R(\xi) = \sinh(\xi/2)/(\xi/2)\) 를 얻었다 (\(R \sim U(-1/2, 1/2]\)). \(R\) 은 \(U - 1/2\) 로 평행이동된 변수이므로 \(\kappa_1(R) = 0\), \(\kappa_r(R) = \kappa_r(U)\) for \(r \geq 2\). 즉 문제 11.3 의 짝수 누율들은 곧 \(K_R\) 의 Taylor 계수다.

\(K_R(\xi) = \log(\sinh(\xi/2)/(\xi/2))\) 를 \(\xi = 0\) 근방에서 전개하면

\[ K_R(\xi) = \frac{\xi^2}{24} - \frac{\xi^4}{2880} + O(\xi^6). \]

누율은 \(K_R(\xi) = \sum_{r \geq 1} \kappa_r \xi^r / r!\) 에서 계수를 읽어낸다.

• \(\xi^2\) 계수 \(= 1/24 = \kappa_2/2!\) → \(\kappa_2 = 1/12\). ✓
• \(\xi^4\) 계수 \(= -1/2880 = \kappa_4/4!\) → \(\kappa_4 = -24/2880 = -1/120\). ✓

두 방법이 일치한다. 이는 11.1 의 주장 2 가 옳다는 간접 증거이기도 하다.

4.3 왜 \(\kappa_4 < 0\) 인가 — 균등분포는 ‘덜 뾰족’

\(\kappa_4\) 는 정규 분포 대비 꼬리의 두꺼움을 잰다 (excess kurtosis). 균등분포는 완전히 평평한 분포 이므로 정규 분포보다 꼬리가 더 얇다 → \(\kappa_4 < 0\). 구체적으로 excess kurtosis 는 \(\kappa_4/\kappa_2^2 = (-1/120)/(1/144) = -6/5 = -1.2\).

문제 11.2 의 \(\epsilon^2/12\) 바닥은 이 \(\kappa_2 = 1/12\) 에서 직접 나온다. 따라서 11.3 은 11.2 의 상수를 구체화 해 주는 필수 보조 정리다.

5 문제 11.4 — Drosophila 발생률의 비선형 적합

5.1 데이터와 모형

Powsner (1935) 는 초파리 Drosophila melanogaster 의 배아(embryonic) 발달 단계 지속 시간을 온도별로 측정했다 (23 배치, Table 8.8). 배치 \(i\) 에서 온도 \(T_i\), 평균 지속 시간 \(y_i\) (시간), 배치 크기 \(m_i\) 다.

McCullagh-Nelder (§8.4.4) 는 유리 함수 (rational function) 온도 모형을 제안한다.

\[ \log \mu = \beta_0 + \beta_1 T + \frac{\beta_{-1}}{T - \delta}, \tag{8.4} \]

• 감마 오차 · 로그 링크 · 가중치 = 배치 크기 \(m_i\).
• \(\delta\) 는 비선형 모수다 — 분모에 \((T - \delta)\) 로 들어가 있어 공변량 자체의 형태를 바꾼다.
• \(\delta\) 가 알려져 있으면 \(u = 1/(T - \delta)\) 를 공변량으로 꽂아 평범한 감마 GLM 이다.

직관: 유리 함수 모형은 무엇을 포착하는가

Arrhenius 법칙 \(\log(\text{rate}) = -\mu/(RT)\) 은 단순 화학 반응 에 대해 성립한다. 그러나 생물학적 발달은 복수의 효소 반응이 직렬로 이어져 있어 Arrhenius 가 예측하는 “로그 속도 vs \(1/T\)” 의 직선이 관측되지 않는다. 데이터는 29-31°C 근방에서 최단 지속 시간을 보이는 U 자형 이다.

유리 함수 \(\beta_{-1}/(T - \delta)\) 는 \(T \to \delta\) 에서 발산하는 점근 을 허용한다. 즉 “어떤 임계 온도 \(\delta\) 에서 효소가 멈춘다” 는 질적 사실을 수식으로 표현한 것이다.

5.2 §11.4 선형화

\(\delta_0\) 근방에서 Taylor 1차 전개:

\[ \frac{1}{T - \delta} \simeq \frac{1}{T - \delta_0} + (\delta - \delta_0) \cdot \frac{1}{(T - \delta_0)^2}. \]

따라서

\[ \log \mu \simeq \beta_0 + \beta_1 T + \beta_{-1} u + \gamma v, \qquad u = \frac{1}{T - \delta_0}, \quad v = \frac{1}{(T - \delta_0)^2}, \quad \gamma = \beta_{-1}(\delta - \delta_0). \]

네 공변량 \((1, T, u, v)\) 로 가중 감마 GLM (log link, weights \(= m_i\)) 을 적합한 뒤 갱신식

\[\delta_1 = \delta_0 + \frac{\hat\gamma}{\hat\beta_{-1}}.\]

출발값은 \(\delta_0 = 0\) (섭씨) 또는 \(\delta_0 = 55\) (시각적 곡선 맞춤) 로 쓴다. 수 회 반복 후 수렴하는 값은 \(\hat\delta \simeq 58.644\) 다.

5.3 수렴 후 결과

네 모수 추정치와 점근 표준오차 (최종 반복에서 \(v \to \beta_{-1} v\) 치환 후):

모수	추정치	SE
\(\beta_0\) (절편)	3.201	1.594
\(\beta_1\) (선형 항)	\(-0.265\)	0.0355
\(\beta_{-1}\) (쌍곡 항)	\(-217.08\)	125.21
\(\delta\) (임계 온도)	58.644	6.48

이탈도는 \(D = 0.32\) on 19 d.f., 잔차 변동 계수 \(\tilde\sigma = (0.32/19)^{1/2} = 0.13\) (13%). Ch.8 에서 이미 본 결과 그대로다.

최소 지속 시간은 \(\partial(\log \mu)/\partial T = 0\) 에서 달성된다.

\[ \beta_1 - \frac{\beta_{-1}}{(T - \delta)^2} = 0 \Rightarrow \hat T = \hat\delta - \left(-\frac{\hat\beta_{-1}}{\hat\beta_1}\right)^{1/2} = 58.64 - (-217.08/-0.265)^{1/2}. \]

계산하면 \((-217.08/-0.265)^{1/2} = (819.2)^{1/2} = 28.62\) 이므로 \(\hat T = 58.64 - 28.62 = 30.01°C\) — 생물학적으로 타당한 값이다.

5.4 프로파일 이탈도 곡선 — 정규 근사의 한계

문제 11.4 의 핵심 요구는 “\(\delta\) 를 50-100°C 범위에서 고정한 채 이탈도를 플롯하여 \(\hat\delta\) 에 대한 정규 근사의 적절성 을 확인하라” 다.

구체적 절차:

1. 격자 \(\delta_k \in \{50, 52, \ldots, 100\}\) 를 잡는다.
2. 각 \(\delta_k\) 에서 \(u_k = 1/(T - \delta_k)\) 를 공변량으로 하는 3 모수 감마 GLM (\(\beta_0, \beta_1, \beta_{-1}\)) 을 적합한다.
3. 이탈도 \(D(\delta_k)\) 를 기록한다.
4. \(D(\delta_k)\) 를 \(\delta_k\) 에 대해 플롯한다.

\(\hat\delta\) 가 정규 근사를 잘 만족하면 \(D(\delta)\) 는 이차 함수 형태여야 한다. 그러나 Fig. 8.8a 는 명확히 비대칭 인 곡선을 보여 준다.

직관: 왜 \(\delta\) 프로파일이 비이차인가

\(\delta\) 는 분모 에 들어간다 — \(1/(T - \delta)\). 이 함수는 \(\delta\) 에 대해 비선형 · 비대칭 이다. \(\delta\) 가 \(\bar T\) 에서 멀어질수록 곡선의 형태가 다르게 변한다.

• \(\delta\) 가 데이터 범위 (14-32°C) 에 가까워지면 적합값이 발산하며 이탈도가 급증한다.
• \(\delta\) 가 멀어지면 (\(\delta \to \infty\)) \(1/(T - \delta) \to 0\) 이 되어 쌍곡 항이 사라지고 모형은 선형 으로 수렴한다. 즉 deviance 가 포화된다.

두 방향이 비대칭이므로 \(\hat\delta\) 의 오차 분포도 비대칭이다. \(\hat\delta = 58.6 \pm 6.48\) 의 Wald 95% 구간은 \([45.9, 71.4]\) 이지만, 프로파일 이탈도 기반 95% 구간은 아래쪽이 더 짧고 위쪽이 더 긴 비대칭 구간이다.

5.5 재모수화 전략

McCullagh-Nelder 는 \(\delta\) 대신 \(\psi = 1/(\delta - \bar T)\) 를 사용하면 이탈도 곡선이 훨씬 이차에 가까워진다고 보고한다 (Fig. 8.8b). 이유는 분명하다.

\[u_i = \frac{1}{T_i - \delta} = \frac{1}{T_i - \bar T - (\delta - \bar T)}.\]

\(\psi = 1/(\delta - \bar T)\) 를 도입하면 이 항을 \(\psi \cdot [\cdots]\) 꼴로 \(\psi\) 에 대해 선형 으로 쓸 수 있다 — 즉 선형 예측자가 \(\psi\) 에 훨씬 덜 뒤틀려 의존한다.

일반 원칙: 공변량의 분모에 들어간 모수 는 원래 스케일에서 비이차 프로파일을 내기 쉽다. 역수 스케일로 재모수화하면 대칭성이 복구된다.

5.6 Python 구현

import numpy as np
import statsmodels.api as sm

# Powsner (1935) Table 8.8 — 23 batches
T = np.array([14.95, 16.16, 16.19, 17.15, 18.20, 19.08, 20.07, 22.14,
              23.27, 24.09, 24.81, 24.84, 25.06, 25.06, 25.80, 26.92,
              27.68, 28.89, 28.96, 29.00, 30.05, 30.80, 32.00])
y = np.array([67.5, 57.1, 56.0, 48.4, 41.2, 37.80, 33.33, 26.50,
              24.24, 22.44, 21.13, 21.05, 20.39, 20.41, 19.45, 18.77,
              17.79, 17.38, 17.26, 17.18, 16.81, 16.97, 18.20])
m = np.array([54, 182, 153, 129, 64, 94, 82, 57, 135, 188,
              217, 141, 37, 84, 196, 104, 148, 83, 95, 232,
              148, 195, 58], dtype=float)

def fit_at_delta(delta):
    u = 1.0 / (T - delta)
    X = np.column_stack([np.ones_like(T), T, u])
    fam = sm.families.Gamma(link=sm.families.links.log())
    res = sm.GLM(y, X, family=fam, freq_weights=m).fit()
    return res

# §11.4 선형화 반복
delta = 55.0
for it in range(15):
    u = 1.0 / (T - delta)
    v = 1.0 / (T - delta) ** 2
    X = np.column_stack([np.ones_like(T), T, u, v])
    fam = sm.families.Gamma(link=sm.families.links.log())
    res = sm.GLM(y, X, family=fam, freq_weights=m).fit()
    beta_m1 = res.params[2]
    gamma   = res.params[3]
    delta_new = delta + gamma / beta_m1
    print(f"iter {it}: delta={delta_new:.4f}, deviance={res.deviance:.4f}")
    if abs(delta_new - delta) < 1e-5:
        break
    delta = delta_new

# 프로파일 이탈도 곡선
delta_grid = np.arange(50, 101, 1.0)
dev_profile = np.array([fit_at_delta(d).deviance for d in delta_grid])

import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 2, figsize=(10, 4))
ax[0].plot(delta_grid, dev_profile, 'o-')
ax[0].axvline(delta, ls='--', color='gray')
ax[0].set_xlabel(r'$\delta$ (°C)'); ax[0].set_ylabel('profile deviance')
ax[0].set_title(r'$\delta$ scale — 비대칭 곡선')

T_bar = T.mean()
psi_grid = 1.0 / (delta_grid - T_bar)
ax[1].plot(psi_grid, dev_profile, 'o-')
ax[1].axvline(1.0/(delta - T_bar), ls='--', color='gray')
ax[1].set_xlabel(r'$\psi = 1/(\delta - \bar T)$'); ax[1].set_ylabel('profile deviance')
ax[1].set_title(r'$\psi$ scale — 더 이차적')
plt.tight_layout(); plt.show()

왼쪽 플롯은 \(\delta\) 스케일에서 비대칭, 오른쪽은 \(\psi\) 스케일에서 거의 이차다. 신뢰구간은 \(\psi\) 스케일에서 Wald 로 계산한 뒤 역변환해서 \(\delta\) 스케일로 옮기면 비대칭 구간이 자연스럽게 나온다.

6 네 문제의 상호 관계

6.1 이론(11.1-11.3) 과 실전(11.4) 의 연결

연결	내용
11.1 → 11.2	반올림 오차의 CGF 구조 → 감마 분산 함수의 구체적 보정
11.3 → 11.2	균등 분포 누율 값 → 11.2 의 \(1/12\) 상수
11.1-11.3 → 11.4	측정 정밀도 가 분산 함수를 바꾼다는 이론 → Table 8.8 의 Duration ± 0.01 기록 정밀도가 모형 적합에 미치는 한계
§11.4 → 11.4	공변량 안 비선형 모수의 선형화 → Powsner 데이터에 \(\delta\) 적합

구체적으로 Table 8.8 의 기록 정밀도를 확인해 보자. 대부분의 지속 시간은 \(\pm 0.01\) 시간 정밀도로 기록되어 있다 — \(\epsilon = 0.01\). 문제 11.2 에 따르면 반올림 바닥은 \(\epsilon^2/12 = 8.33 \times 10^{-6}\). 배치 평균의 분산은 \((\mu^2/\nu)/m_i\) 수준으로, \(\mu = 20\) · \(m_i = 100\) · \(\tilde\sigma = 0.13\) 을 대입하면 \(0.169/100 = 1.69 \times 10^{-3}\). 반올림 바닥은 이보다 200 배 작으므로 무시 가능하다. 감마 모형이 정당화된다.

6.2 일반 원칙 — “이상 모형” 과 “기록된 모형” 의 차이

11.1 의 Kolassa-McCullagh 정리가 주는 전체적 교훈은 다음과 같다.

“기록된 값 \(Y\) 는 이상적 관측값 \(Z\) 와 누율 수준에서 홀수 누율만 같고, 짝수 누율은 \(\epsilon^r \kappa_r(R)\) 만큼 더해진다. 이 차이가 추론에 영향을 미치는 범위는 \(\epsilon^r\) 에 비례해 빠르게 작아진다.”

이는 “유한 정밀도 기록이 통계적 모형을 어떻게 왜곡시키는가” 에 대한 정량적 경계 를 제공한다. 실무에서는 이 경계를 계산해 보고 나서야 “데이터 정밀도가 추정치에 미치는 영향” 을 감당할 수준인지 판단할 수 있다.

7 관련 주제

선행 지식

• Non-Linear Parameters in the Covariates — Box-Tidwell 선형화 (McCullagh §11.4)
• Non-Linear Parameters — Examples (McCullagh §11.5)
• The Gamma Distribution — 밀도·누율·형태 모수 (McCullagh §8.2)
• Gamma GLM Examples — 자동차 보험·혈액 응고·강수량·초파리 발생률 (McCullagh §8.4)

직접 관련 — Ch.11 시리즈 다른 연습 장

• Ch.10 Exercises — 지수 이탈도·공통 모수·Table 10.1 유도 (McCullagh §10.9)
• Quasi-likelihood 심화 결과와 연습 — Poisson-Gamma 혼합·Löwner 순서 (McCullagh §9.8)

관련 개념

• 모멘트 생성 함수와 누율 (Casella-Berger Ch.2.3)
• 균등 분포 — 연속형 (McCullagh §3 및 Casella-Berger)
• 점근 정규성과 프로파일 우도

후속 주제

• Model Checking — 스코어 검정·레버리지·영향력 (McCullagh Ch.12) — 비선형 모수 추정 후의 진단 전략