개인적으로 Profile Likelihood(PL)를 공부하다가 PL을 이용한 신뢰구간을 어떻게 계산하는지 궁금해서 찾아본 논문이 'A Method for Computing Profile-Likelihood-Based Confidence Intervals'이었다.
이번 포스팅에서는 이 논문을 읽고 요약한 내용을 정리하려고 한다.
1. 'Profile-Likelihood-Based' Confidence Intervals
본 포스팅에서는 수식을 포함하고 있습니다.
티스토리 피드에서는 수식이 제대로 표시되지 않을 수 있으니
PC 웹 브라우저 또는 모바일 웹 브라우저에서 보시기 바랍니다.
1. 'Profile-Likelihood-Based' Confidence Intervals
전통적으로 신뢰구간을 계산하는 방법은 최대 우도 추정량 $\hat {\theta}$의 점근 정규분포 성질의 기초하고 있다. 하지만 최대 우도 추정량은 소표본인 경우 그 분포가 정규분포와 매우 달라지는 등, 그 성질이 매우 달라질 수 있다. 좀 더 강건한(Robust) 신뢰구간 계산방법으로 우도비의 점근 카이제곱 분포 성질을 이용하는 방법이 있다.
$l(\theta)$를 $\theta \subset R^k$의 로그 우도 함수라고 하자. 이때 점근 $theta$의 $1-\alpha$ 신뢰 영역은 다음과 같다.
$$\{ \theta : 2[l(\hat{\theta}) - l(\theta)] \leq q_k(1-\alpha) \}$$
여기서 $q_k(1-\alpha)$은 자유도가 $k$인 카이제곱 분포에서 $1-\alpha$ 분위수이다.
이때 $\theta$의 원소($\theta$는 벡터이다) 중에서 관심에 대상이 되는 모수를 $\beta=\theta_j$, $\Theta_j(\beta) = \{ \theta : \theta_j=\beta \}$ 라 하자.
$$l_j(\beta) = \max_{\theta \in \Theta_j(\beta)} \tag{1}$$
을 $\beta$에 대한 Profile 우도 함수라고 한다.
이때 $\beta$에 대한 점근 $1-\alpha$ Profile-Likelihood-Based(PLB) 신뢰구간은 다음과 같다.
$$\{ \theta : 2[l(\hat{\theta}) - l_j(\beta)] \leq q_1(1-\alpha) \}$$
$w = (\theta_1, \ldots, \theta_{j-1}, \theta_{j+1}, \ldots, \theta_k)^t$라 하자. 즉, $w$는 관심의 대상이 되는 모수를 제외한 나머지 모수 벡터를 의미한다. 이때 $\hat{w}(\beta)$를 (1)의 최대값이 되도록 하는 벡터이며 이는 $\Theta_j(\beta)$의 Interior에 존재한다고 가정한다.
지금까지는 PLB 신뢰구간을 효율적으로 계산하는 방법에 대해서 무관심했다. 왜냐하면 기존에 모수 하나씩 신뢰구간을 구하는 방법이 있었기 때문이다. 하지만 $k$차원 모수에 대해서 하나씩 계산하면 총 $2k$의 계산량이 들어가게 된다(하나의 모수에 대하여 신뢰 하한, 상한 두 개를 구해야 하므로).
하지만 이 논문에서는 다음의 $k$개 방정식을 풀면 신뢰구간의 양 끝 값을 한 번에 구할 수 있다고 한다.
$$\begin{equation*}
\begin{pmatrix}
l(\theta) - l^* \\
\frac{\partial l}{\partial w}(\theta)
\end{pmatrix}=0
\end{equation*}\tag{2}$$
이 시스템의 해는 각 스텝에서 $\hat{w}(\beta)$상의 벡터의 중간값을 필요로 하지 않는다고 한다(? 뭔 말인지는 읽어보면 알 수 있겠지?). 이러한 점이 기존 alternative 하게 업데이트하는 과정보다 계산상의 이점이 있다고 한다. PLB의 신뢰구간 양 끝 값 계산이 원래 로그 우도 함수의 최대 우도 추정량을 구하는 것과 비슷한 난이도라고 한다(이러면 PLB 쓰는 이유가 없는 거 아닌가?).
2. Basis of the Algorithm
먼저 이 논문에서는 $\beta=\theta_1$으로 고정시킨다. 그리고 $\beta, w$를 초기화하거나, 로그 우도 함수 $l(\theta)$의 이차 근사 계산을 할 때에는 $\hat{\theta}$는 주어졌다고 가정한다.
$\hat{\theta}$로부터 바로 $(\beta_L, \hat{w}(\beta_L))$ 또는 $(\beta_U, \hat{w}(\beta_U))$을 유도한다면 그 정확성을 충분하지 않을 수 있다고 한다. $\beta_L, \beta_U$는 각각 $\beta$ 신뢰구간의 하한과 상한이다. (2)의 해를 찾기 위해 기본적인 Newton-Raphson 알고리즘은 다음과 같다.
$$\begin{equation*}
\begin{pmatrix}
\beta^{(i+1)}\\
w^{(i+1)}
\end{pmatrix}=\begin{pmatrix} \beta^{(i)} \\ w^{(i)} \end{pmatrix} - \begin{pmatrix} \frac{\partial l}{\partial \beta} & \left (\frac{\partial l}{\partial w} \right )^t \\ \frac{\partial^2 l}{\partial\beta \partial w} & \frac{\partial^2 l}{\partial w^2} \end{pmatrix}^{-1}\begin{pmatrix} l-l^* \\ \frac{\partial l}{\partial w} \end{pmatrix} \end{equation*}\tag{3}$$
이때 (3) 안에 행렬 부분은 첫 번째 행이 0이어서 singular하다. 이를 해결하기 위해 몇 가지 고려사항이 필요하다. 기존 Newton-Raphson 알고리즘을 수정하여 해를 구하는 것이 일반적으로 이득이 된다고 한다.
수정 방법은 다음과 같다. 다음과 같이 스텝 사이즈 $h$를 도입한다.
$$\begin{equation*} \theta^{(1)} = \begin{pmatrix} \beta^{(1)} \\ w^{(1)} \end{pmatrix} = \begin{pmatrix} \hat{\beta} \\ \hat{w}(\hat{\beta}) \end{pmatrix} + h \begin{pmatrix} 1 \\ \frac{d \hat{w}}{d\beta} \end{pmatrix}\end{equation*}$$
이때
$$d\hat{w}/d\beta = -\left (\frac{\partial^2 l}{\partial w^2} \right )^{-1}\frac{\partial^2 l}{\partial \beta \partial w}$$
이다.
여기서 $h$의 선택에 대한 문제가 생기는데 여러 기준을 이용하여 h를 선택할 수 있다고 한다.
저자는 수정된 Newton-Raphson 방법을 제안한다. 이 방법은 $l$의 이차 미분을 고려한다고 한다. 이렇게 수정하면 수렴 속도가 더 빠르다. 그리고 수정된 알고리즘을 통하여 로그 함수 평면이 실제로 이차식으로 잘 근사가 되는지를 확인할 수 있다고 한다. 이때 잘 근사가 안된다면 수정 버전이 아닌 기존 Newton-Raphson 버전을 이용하면 된다고 한다.
3. Performance in Simulations
만약 $l(\theta)$가 이차 형태라면 단 두 스텝만에 신뢰구간을 계산할 수 있다고 한다. 여기서는 두 가지 유형의 로그 우도 함수에 대하여 실험했다고 한다.
첫 번째 경우 큐빅 형태의 로그 우도 함수에 대해서 실험했는데 평균적으로 6.3 스텝 만에 수렴했다고 한다.
와이블 분포에 대해서도 평균적으로 8 스텝만에 수렴했다고 한다.
4. Examples
실제 두 데이터에 대해서 저자가 소개한 알고리즘을 적용하여 신뢰구간을 계산한다.
특히 두 번째 데이터의 경우 기존 Newton-Raphson 알고리즘으로는 신뢰구간을 찾을 수 없으며 논문에서 제안한 수정된 Newton-Raphson 방법을 이용하면 안정적으로 신뢰구간을 찾을 수 있음을 밝혔다.
5. Discussion
수정된 Newton-Raphson 알고리즘의 수렴성은 로그 우도 함수가 얼마나 이차식과 잘 근사 하느냐에 의존한다. 이 수정된 Newton-Raphson 알고리즘은 초기값을 잘 선정해야 하는데 정교한 방법을 이용하면 그게 가능하다고 한다. 로그 우도 방정식이 이차식과 매우 다른 경우 스텝 사이즈 $h$를 줄이거나 미분방정식 테크닉을 이용하면 좀 더 나은 초기값을 찾을 수 있다고 한다.
식 (2)의 솔루션을 찾는데 Newton-Raphson 이외의 방법을 적용할 수 있다고 한다. 관측된 $l$의 이차 미분이 아닌 기대값을 적용할 수도 있다고 한다.
또한 계산량이 많이 필요한 붓스트랩이나 분산 안정 변환 등의 방법을 이용하여 신뢰구간을 계산할 수도 있다고 한다. 그리고 최대 우도 추정량에 대하여 신뢰구간이 대칭성을 만족하도록 하는 모수 변환 방법도 있으며 이러한 변환 방법을 이용하면 더 빠르게 (변환된 모수에 대한) 신뢰구간을 더 빠르게 얻을 수 있다고 한다.
댓글