Models with Gamma-Distributed Observations

분산 함수·이탈도·정준/로그/항등 연결·산포 추정의 완전한 전개 (McCullagh §8.3)

감마 GLM의 모형 구성 요소를 상세히 다룬다. 분산 함수 \(V(\mu)=\mu^2\) 의 지수족 유도, 이탈도와 \(D^+\) 대안, 정준(역수) 연결의 역다항식 반응곡면, 로그 연결의 곱셈적 모형과 가중 함수 동일성, 항등 연결의 분산 성분 추정, 산포 모수의 MLE(디감마 방정식)와 모멘트 추정량까지 McCullagh §8.3 전체를 전개한다.

Statistics
GLM
저자

Kwangmin Kim

공개

2026년 04월 18일

1 개관 — §8.3의 위치

§8.2에서 감마 분포 \(G(\mu,\nu)\) 의 밀도, 누율, 합성 성질을 다루었다 — 특히 변동계수 일정 조건 \(\operatorname{var}(Y) = \sigma^2 \mu^2\), 즉 분산 함수 \(V(\mu) = \mu^2\) 이 감마를 양의 연속 자료의 기본 분포로 선택하게 한 핵심 근거였다. 이 절에서는 그 분산 함수를 GLM 프레임워크에 올려놓는다. 구체적으로:

  • 지수족 표현에서 분산 함수와 정준 모수를 도출하고 (§8.3.1)
  • 이탈도를 유도하며 (§8.3.2)
  • 세 가지 연결 함수(역수·로그·항등)에 따른 모형 계열을 상세히 비교하고 (§8.3.3-§8.3.5)
  • 산포 모수 \(\sigma^2 = 1/\nu\) 의 추정법을 제시한다 (§8.3.6)

직관. Ch.3(정규), Ch.4(이항), Ch.6(포아송)에서 보았듯 GLM은 항상 세 구성 요소 – 확률 성분(분산 함수), 체계적 성분(선형 예측자), 연결 함수 – 로 구성된다. 감마 GLM도 같은 틀이지만, \(V(\mu)=\mu^2\) 라는 이차 분산 함수가 역수·로그·항등 연결 각각에서 질적으로 다른 모형 계열을 만든다는 점이 핵심이다.

2 분산 함수

2.1 지수족 표준형으로부터의 유도

감마 로그 우도를 \(\mu\) 와 \(\nu\) 의 함수로 쓰면

\[ \ell(\mu,\nu;y) = \nu\!\left(-\frac{y}{\mu} - \log\mu\right) + \nu\log y + \nu\log\nu - \log\Gamma(\nu). \]

Ch.2의 지수족 표준형 \(\ell = \{y\theta - b(\theta)\}/a(\phi) + c(y,\phi)\) 과 대조하면

지수족 요소 감마 대응
정준 모수 \(\theta\) \(-1/\mu\)
누율 함수 \(b(\theta)\) \(-\log(-\theta) = \log\mu\)
산포 함수 \(a(\phi)\) \(1/\nu = \sigma^2\)

직관. \(\theta = -1/\mu\) 라는 것은, 감마 분포의 “자연스러운 좌표”가 평균의 음의 역수라는 뜻이다. 이것이 정준 연결 \(\eta = 1/\mu\) (부호 규약에 따라 \(-1/\mu\))를 결정한다.

\(b(\theta)\) 를 미분하면

\[ b'(\theta) = \frac{-1}{\theta} = \mu, \qquad b''(\theta) = \frac{1}{\theta^2} = \mu^2. \]

따라서 분산 함수는

\[ \boxed{V(\mu) = b''(\theta) = \mu^2.} \]

직관. 분산 함수 \(V(\mu) = \mu^2\) 는 “절대 오차의 크기가 평균에 비례한다”는 일상 경험을 수학적으로 포착한다. 10만 원짜리 보험금에서 1만 원의 오차와, 100만 원짜리 보험금에서 1만 원의 오차는 의미가 다르다 – 상대적으로 보면 전자가 10배 더 심각하다. \(V(\mu) = \mu^2\) 는 이 상대적 관점을 자동으로 반영한다.

2.2 GLM 분산 함수 계층

분포 \(V(\mu)\) 멱지수 \(p\) 해석
정규 \(1\) \(0\) 절대 오차 일정
포아송 \(\mu\) \(1\) SD \(\propto \sqrt{\mu}\)
감마 \(\mu^2\) \(2\) CV 일정
역가우스 \(\mu^3\) \(3\) SD \(\propto \mu^{3/2}\)

Tweedie 족 \(V(\mu) = \mu^p\) 에서 감마는 \(p = 2\) 에 해당한다. \(p\) 가 커질수록 분산의 평균 의존도가 강해진다.

3 이탈도(Deviance)

3.1 유도

\(\nu\) 를 알려진 상수로 취급하면, 독립 관측의 로그 우도는

\[ \ell(\mu) = \sum_i \nu\!\left(-\frac{y_i}{\mu_i} - \log\mu_i\right). \]

가중 관측(\(\nu_i = \nu w_i\), \(w_i\) 는 알려진 가중치)인 경우

\[ \ell(\mu) = \nu\sum_i w_i\!\left(-\frac{y_i}{\mu_i} - \log\mu_i\right). \]

포화 모형(\(\hat\mu_i = y_i\))의 로그 우도는

\[ \ell_{\max} = -\nu\sum_i w_i(1 + \log y_i), \]

\(y_i > 0\) 일 때만 유한하다.

이탈도 \(D = 2[\ell_{\max} - \ell(\hat\mu)]/\nu\) 를 정리하면

\[ \boxed{ D(y;\hat\mu) = -2\sum_i w_i \left\{\log\!\left(\frac{y_i}{\hat\mu_i}\right) - \frac{y_i - \hat\mu_i}{\hat\mu_i}\right\}. } \]

직관. 이탈도의 두 항을 분리해서 읽으면:

  • \(\log(y_i/\hat\mu_i)\): 비율 편차. 관측이 적합의 2배이면 \(\log 2 \approx 0.69\)
  • \((y_i - \hat\mu_i)/\hat\mu_i\): 상대 잔차. 관측이 적합의 2배이면 값은 1

정규 모형의 이탈도 \(\sum(y_i - \hat\mu_i)^2\) 가 절대 편차의 제곱을 재는 것과 달리, 감마 이탈도는 비율 편차를 잰다. 이는 변동계수가 일정한 자료의 본질에 정확히 부합한다.

3.2 절편이 있는 모형에서의 단순화

모형에 절편이 포함되면 \(\sum w_i y_i/\hat\mu_i = \sum w_i\) 이므로

\[ \sum_i w_i \frac{y_i - \hat\mu_i}{\hat\mu_i} = 0. \]

따라서 이탈도의 마지막 항이 사라지고

\[ D(y;\hat\mu) = -2\sum_i w_i \log\!\left(\frac{y_i}{\hat\mu_i}\right) \]

만 남는다 (Nelder & Wedderburn, 1972).

직관. 절편이 있으면 적합값의 가중 평균이 관측값의 가중 평균과 일치한다. 이 조건이 상대 잔차의 합을 0으로 만들어, 이탈도가 순수하게 비율 편차만의 함수가 된다.

3.3 \(y_i = 0\) 인 경우: 대안 통계량 \(D^+\)

\(y_i = 0\) 이면 \(\log(y_i/\hat\mu_i) = -\infty\) 로 이탈도가 정의되지 않는다. 이 경우 대안 통계량

\[ D^+(y;\hat\mu) = 2C(y) + 2\sum_i w_i \log\hat\mu_i + 2\sum_i w_i y_i/\hat\mu_i \]

를 사용한다. \(C(y)\) 는 \(y\) 만의 임의 유계 함수이다.

\(D\) 와 \(D^+\) 의 비교:

\(D(y;\hat\mu)\) \(D^+(y;\hat\mu)\)
양수성 항상 \(\ge 0\) 보장 안 됨
\(y_i = 0\) 허용 불가 (\(-\infty\)) 가능
\(\hat\nu\) 의 MLE 함수 관계 있음 없음

직관. \(D^+\) 는 \(\log y_i\) 항을 \(C(y)\) 속에 흡수하여 영의 관측을 허용한다. 대가로 잔차 제곱합 같은 직관적 해석과, \(\nu\) 의 MLE와의 연결을 잃는다. 실무에서 반올림 오차로 인한 가짜 영(spurious zero)이 문제가 되는 경우에 유용하다.

7 산포 모수의 추정

감마 GLM에서 모수 추정량 \(\hat\beta\) 의 근사 공분산 행렬은

\[ \operatorname{cov}(\hat\beta) \simeq \sigma^2(X^TWX)^{-1}, \]

\[ W = \operatorname{diag}\!\left\{ \frac{(d\mu_i/d\eta_i)^2}{V(\mu_i)} \right\} \]

이므로 \(\sigma^2\) (또는 \(\nu = 1/\sigma^2\))의 추정이 필수적이다.

7.1 MLE: 디감마 방정식

감마 모형 하에서 \(\nu\) 의 MLE는 다음 방정식의 해이다.

\[ \boxed{ 2n\{\log\hat\nu - \psi(\hat\nu)\} = D(y;\hat\mu), } \tag{8.1} \]

여기서 \(\psi(\nu) = \Gamma'(\nu)/\Gamma(\nu)\) 는 디감마 함수이다.

유도의 핵심. 프로파일 로그 우도 \(\ell(\nu) = \sum[\nu\log\nu - \log\Gamma(\nu) + (\nu-1)\log y_i - \nu y_i/\hat\mu_i - \nu\log\hat\mu_i]\) 를 \(\nu\) 에 대해 미분하면 \(\log\nu - \psi(\nu) + \bar{\log y} - \overline{y/\hat\mu} - \overline{\log\hat\mu} = 0\) 이 되고, 이탈도의 정의를 대입하면 식 (8.1)이 된다.

직관. 좌변의 \(g(\nu) = \log\nu - \psi(\nu)\) 는 \(\nu\) 의 단조 감소 함수이므로 해가 유일하다.

  • \(g(\nu) \to \infty\) (\(\nu \to 0^+\))
  • \(g(\nu) \approx 1/(2\nu)\) (\(\nu \to \infty\))

따라서 \(\hat\nu\) 가 클 때(작은 \(\sigma^2\)) \(D \approx n/\hat\nu = n\sigma^2\) 으로, 정규 모형에서 “잔차 제곱합 \(\approx n\sigma^2\)”인 관계의 감마 버전이다.

7.2 편향 보정

\(p\) 개 모수를 추정한 효과를 반영하면

\[ 2n\{\log\hat\nu - \psi(\hat\nu)\} - p\hat\nu^{-1} = D(y;\hat\mu). \tag{8.2} \]

보정항 \(-p\hat\nu^{-1}\) 은 \(E[D(Y;\hat\mu)]\) 의 점근 전개에서 \(O(1)\) 항이다. 정규 이론에서 \(\hat\sigma^2 = \text{RSS}/n\) (MLE)과 \(s^2 = \text{RSS}/(n-p)\) (비편향)의 관계와 명확한 유추가 된다.

7.3 큰 \(\nu\) 에 대한 근사

\(\nu\) 가 충분히 크면 \(g(\nu) \approx 1/(2\nu) + 1/(12\nu^2) + \cdots\) 이므로 \(O(\nu^{-2})\) 이상을 무시하면

\[ \hat\nu^{-1} \simeq \frac{\bar{D}(6 + \bar{D})}{6 + 2\bar{D}}, \qquad \bar{D} = D(y;\hat\mu)/n. \]

직관. 이 근사는 디감마 함수의 수치 풀이 없이 이탈도만으로 \(\sigma^2\) 를 빠르게 추정할 수 있게 한다. \(\bar{D}\) 가 작으면(좋은 적합) \(\hat\nu^{-1} \approx \bar{D}\) 로, “평균 이탈도 \(\approx\) 산포”라는 단순한 관계가 성립한다.

7.4 MLE의 문제점

이탈도 기반 추정량의 근본적 문제 두 가지:

  • 영 민감도: \(y_i = 0\) 이면 \(D = \infty\), \(\hat\nu = 0\). 반올림 오차로 인한 가짜 영에도 취약하다.
  • 비일치성: 감마 가정이 틀리면 \(\hat\nu^{-1}\) 은 변동계수를 일치 추정하지 못한다. 분포 형태가 감마에서 벗어나면 \(\log(y/\hat\mu)\) 항의 기댓값이 \(-(y-\hat\mu)/\hat\mu\) 와 맞지 않기 때문이다.

7.5 모멘트 추정량: Pearson \(X^2\) 기반

McCullagh & Nelder가 선호하는 대안은

\[ \boxed{ \tilde\sigma^2 = \frac{1}{n-p}\sum_i \left(\frac{y_i - \hat\mu_i}{\hat\mu_i}\right)^{\!2} = \frac{X^2}{n-p}. } \tag{8.3} \]

여기서 \(X^2\) 는 Pearson 카이제곱 통계량이다.

직관. 이 추정량은 “상대 잔차의 제곱 평균”이다. 정규 이론의 \(s^2 = \sum(y_i - \hat\mu_i)^2/(n-p)\) 가 절대 잔차를 쓰는 것에 대응하여, 감마의 \(\tilde\sigma^2\) 는 상대 잔차를 쓴다. 변동계수가 일정한 자료에서 절대 잔차보다 상대 잔차가 더 자연스러운 척도이므로, 이 선택은 자료의 본성에 부합한다.

핵심 장점:

  • 감마 가정 없이도 \(\sigma^2\) 를 일치 추정한다 – \(\beta\) 의 일치 추정만 있으면 충분하다
  • \(y_i = 0\) 이어도 정의된다
  • 계산이 간단하다

7.6 모멘트 추정량의 편향

단일 감마 표본에서

\[ E(\tilde\sigma^2) = \sigma^2\!\left[1 - \frac{\sigma^2}{n} + O(n^{-2})\right]. \]

음의 편향이 존재하며, 제수를 \(n - p\) 로 해도 \(O(n^{-1})\) 편향이 완전히 제거되지 않는다.

직관. 정규 이론에서 \(s^2\) 가 비편향인 것은 \(V(\mu) = 1\) (상수)이기 때문이다. \(V''(\mu) > 0\) (볼록한 분산 함수)이면 Jensen 부등식에 의해 \(E[V(\hat\mu)] > V(E[\hat\mu])\) 가 되어 분모가 과대평가되고, 결과적으로 \(\tilde\sigma^2\) 가 과소추정된다. 감마의 \(V(\mu) = \mu^2\) 은 볼록이므로 음의 편향이 발생하는 것이다.

7.7 두 추정량의 비교 요약

MLE (\(\hat\nu^{-1}\)) 모멘트 (\(\tilde\sigma^2\))
감마 가정 필요 필수 불필요
\(y_i = 0\) 허용 불가 가능
편향 \(O(n^{-1})\), 보정 가능 (식 8.2) \(O(n^{-1})\), 보정 불완전
효율 (감마 참) 최적 약간 낮음
견고성 (감마 위반) 비일치 일치

8 Python 코드 예시

8.1 Step 1: 감마 GLM 직접 구현 — 역수 연결

import numpy as np

def gamma_glm_irls(X, y, w=None, link='reciprocal', max_iter=25, tol=1e-8):
    """감마 GLM의 IRLS 알고리즘 (순수 numpy 구현)."""
    n, p = X.shape
    if w is None:
        w = np.ones(n)

    # 초기값: mu = y (포화 모형에서 출발)
    mu = y.copy().astype(float)
    mu[mu <= 0] = 0.1  # 영 방지

    for iteration in range(max_iter):
        # 연결 함수에 따른 eta, d_mu/d_eta
        if link == 'reciprocal':
            eta = 1.0 / mu
            dmu_deta = -mu**2        # d(1/eta)/deta = -1/eta^2
        elif link == 'log':
            eta = np.log(mu)
            dmu_deta = mu
        else:  # identity
            eta = mu.copy()
            dmu_deta = np.ones(n)

        # V(mu) = mu^2 (감마 분산 함수)
        V = mu**2

        # IRLS 가중치: W = w * (dmu/deta)^2 / V(mu)
        W = w * dmu_deta**2 / V

        # 조정 종속 변수: z = eta + (y - mu) / (dmu/deta)
        z = eta + (y - mu) / dmu_deta

        # 가중 최소제곱: (X'WX) beta = X'Wz
        XtWX = X.T @ np.diag(W) @ X
        XtWz = X.T @ (W * z)
        beta = np.linalg.solve(XtWX, XtWz)

        # eta, mu 업데이트
        eta_new = X @ beta
        if link == 'reciprocal':
            mu_new = 1.0 / eta_new
        elif link == 'log':
            mu_new = np.exp(eta_new)
        else:
            mu_new = eta_new.copy()

        # 수렴 판정
        if np.max(np.abs(mu_new - mu) / (np.abs(mu) + 1e-10)) < tol:
            mu = mu_new
            break
        mu = mu_new

    # 이탈도
    dev = -2 * np.sum(w * (np.log(y / mu) - (y - mu) / mu))

    # 모멘트 추정량
    sigma2_moment = np.sum(w * ((y - mu) / mu)**2) / (n - p)

    return beta, mu, dev, sigma2_moment


# 예시: 혈액 응고 데이터 (Lot 1)
u = np.array([5, 10, 15, 20, 30, 40, 60, 80, 100], dtype=float)
y = np.array([118, 58, 42, 35, 27, 25, 21, 19, 18], dtype=float)
x = np.log(u)  # 로그 스케일 공변량

X = np.column_stack([np.ones(len(x)), x])
beta, mu_hat, dev, sigma2 = gamma_glm_irls(X, y, link='reciprocal')

print(f"beta = {beta}")
print(f"Deviance = {dev:.4f}")
print(f"sigma^2 (moment) = {sigma2:.4f}")
print(f"CV = {np.sqrt(sigma2):.4f}")

8.2 Step 2: statsmodels 활용

import statsmodels.api as sm
import numpy as np

u = np.array([5, 10, 15, 20, 30, 40, 60, 80, 100], dtype=float)
y = np.array([118, 58, 42, 35, 27, 25, 21, 19, 18], dtype=float)
x = np.log(u)
X = sm.add_constant(x)

# 역수 연결 + 감마
model_inv = sm.GLM(y, X, family=sm.families.Gamma(sm.families.links.InversePower()))
result_inv = model_inv.fit()
print("=== 역수 연결 ===")
print(result_inv.summary())
print(f"Deviance = {result_inv.deviance:.4f}")
print(f"Pearson sigma^2 = {result_inv.pearson_chi2 / result_inv.df_resid:.4f}")

# 로그 연결 + 감마
model_log = sm.GLM(y, X, family=sm.families.Gamma(sm.families.links.Log()))
result_log = model_log.fit()
print("\n=== 로그 연결 ===")
print(result_log.summary())
print(f"Deviance = {result_log.deviance:.4f}")
print(f"Pearson sigma^2 = {result_log.pearson_chi2 / result_log.df_resid:.4f}")

8.3 R 코드 참고

# 혈액 응고 데이터 (Lot 1)
u <- c(5, 10, 15, 20, 30, 40, 60, 80, 100)
y <- c(118, 58, 42, 35, 27, 25, 21, 19, 18)
x <- log(u)

# 역수 연결 (정준)
fit_inv <- glm(y ~ x, family = Gamma(link = "inverse"))
summary(fit_inv)

# 로그 연결
fit_log <- glm(y ~ x, family = Gamma(link = "log"))
summary(fit_log)

# 산포 추정: 모멘트 vs MLE
cat("Pearson sigma^2 =", summary(fit_inv)$dispersion, "\n")
cat("Deviance/(n-p)  =", fit_inv$deviance / fit_inv$df.residual, "\n")

9 세 연결 함수의 비교 요약

역수 \(\eta = 1/\mu\) 로그 \(\eta = \log\mu\) 항등 \(\eta = \mu\)
지수족 지위 정준 (충분통계량 선형) 비정준 비정준
모형 해석 속도(rate), 역다항식 곱셈적(multiplicative) 분산 성분
\(\eta\) 범위 \((0, \infty)\) – 양수 제약 \((-\infty, \infty)\) – 제약 없음 \((0, \infty)\) – 양수 제약
IRLS 가중치 \(\mu^2/\mu^2 \cdot (1/\mu^4) \cdot \mu^4\) \(= 1\) (부분적) \(\mu^2/\mu^2 = 1\) (정확히) \(1/\mu^2\)
전형적 응용 생물학 반응곡면, 보험 경제학, 임상 시간 ANOVA 분산 성분
수치 안정성 주의 필요 (\(\hat\mu < 0\) 위험) 우수 주의 필요 (\(\hat\mu < 0\) 위험)

실무 선택 지침. 특별한 이유가 없다면 로그 연결이 기본 선택이다. 양수 제약이 자동으로 충족되고, 가중치가 1이라 수치적으로 안정적이며, 모수가 곱셈적 효과로 직관적으로 해석된다. 역수 연결은 반응곡면 모형에서, 항등 연결은 분산 성분 추정에서 각각 자연스러운 선택이다.

Subscribe

Enjoy this blog? Get notified of new posts by email: