본 포스트는 필자가 추정 이론을 공부한 내용을 바탕으로 정리한 포스팅이다.
1. Introduction
추정 이론(estimation theory)는 관측된 데이터를 바탕으로 모델의 파라미터나 상태를 예측하는 다양한 방법을 정리한 이론이다. 이는 데이터 분석, 신호처리, 기계학습, 금융, 로봇공학 등 다양한 분야에서 널리 쓰이고 있으며 주로 불확실성을 다루는 과정에서 정확한 결정을 내리기 위한 필수적인 도구로 사용된다. 추정 이론의 응용 분야는 매우 넓은데 통신에서는 신호의 품질을 추정하거나 기계학습에서는 데이터로부터 알고리즘의 파라미터를 결정하는데 사용된다. 또한 금융 분야에서는 시장의 미래 동향을 예측하기 위한 변수를 추정하는데 필수적으로 사용되고 있다.
1.1. The Mathematical Estimation Problem
다음과 같은 데이터 모델이 주어졌다고 가정하자.
\begin{equation}
\begin{aligned}
x[n] = \theta + w[n] \quad n=0,1,\cdots,N-1
\end{aligned}
\end{equation}
- $x[n]$: 관측된 데이터
- $\theta$: 추정하고자 하는 미지의 파라미터
위 식에서 $w[n]$은 랜덤 노이즈라고 하며 일반적으로 평균이 $0$이고 분산은 $\sigma^2$인 가우시안 분포를 따른다고 가정한다.
\begin{equation}
\begin{aligned}
w[n] \sim \mathcal{N}(0, \sigma^{2})
\end{aligned}
\end{equation}
- $a \sim \mathcal{N}(\mu, \sigma^{2})$ : 확률변수 $a$가 평균이 $\mu$이고 분산이 $\sigma^{2}$인 가우시안 분포를 따른다.
좋은 추정값(estimator)을 얻기 위해서는 우선 수학적으로 데이터를 잘 모델링해야 한다. 데이터는 랜덤성을 띄기 때문에 확률 밀도 함수(probability density function, pdf) $p(x[0], x[1], \cdots, x[N-1]; \theta)$를 사용하여 데이터를 표현한다. 이 때 $x[n]$은 $N$개의 데이터를 의미하며 $\theta$는 미지의 모델 파라미터를 의미한다. 만약 $N=1$인 경우 pdf는 아래와 같이 나타낼 수 있다.
\begin{equation}
\boxed{ \begin{aligned}
p(x[0]; \theta) = \frac{1}{ \sqrt{ 2\pi\sigma^{2} }} \exp\bigg[ -\frac{1}{2\sigma^{2}} (x[0] - \theta)^{2} \bigg]
\end{aligned} }
\end{equation}
pdf를 표기하는 서로 다른 방법을 정리하면 다음과 같다.
1. $p(x; \theta)$ : 매개 변수 $\theta$에 대한 $x$의 확률 분포 (pdf of $x$ paramterized by $\theta$)
2. $p(x|\theta)$ : $\theta$가 주어졌을 때 $x$의 조건부 확률 분포 (pdf of $x$ given $\theta$)
3. $p(x,\theta)$ : $x, \theta$가 동시에 발생할 결합 확률 분포 (joint pdf of $x, \theta$)
1.에서 $\theta$는 확률 변수(random variable)일수도 있고 아닐수도 있다. $x$는 확률 변수이다.
2.,3.에서 $x$와 $\theta$는 반드시 확률 변수이어야 한다.
3. 에서 두 확률 변수 $x, \theta$는 일반적으로 서로 독립이 아니며 종속적이다.
위 식에서 보다시피 파라미터 $\theta$는 $x[0]$의 확률에 직접적인 영향을 주기 때문에, 역으로 $x[0]$의 값을 보고 $\theta$를 추정하는 것이 가능하다.
예를 들면 아래와 같은 그림이 주어졌다고 했을 때 $x[0]$ 값이 음수가 관측되면 우리는 $\theta=\theta_{2}$보다는 $\theta=\theta_1$이라고 일반적으로 유추할 수 있다. 하지만 실제 문제에서는 위와 같은 pdf는 주어지지 않기 때문에 데이터 $x$와 파라미터 $\theta$의 관계를 정의해야 한다.
다음으로 다른 유형의 데이터 모델을 살펴보자.
\begin{equation}
\begin{aligned}
x[n] = A + Bn + w[n] \quad \quad n=0,1,\cdots,N-1
\end{aligned}
\end{equation}
이 때, 추정하고자 하는 파라미터는 $\theta = [A,B]^{\intercal}$가 된다. $N$개의 데이터를 벡터화하여 $\mathbf{x} = [x[0], x[1], \cdots, x[N-1]]^{\intercal}$라고 표현하면 $\mathbf{x}$에 대한 pdf는 다음과 같다.
\begin{equation} \label{eq:intro1}
\boxed{ \begin{aligned}
p(\mathbf{x}; \theta) = \frac{1}{ (2\pi\sigma^{2})^{\frac{N}{2}} } \exp\bigg[ -\frac{1}{2\sigma^{2}} \sum_{n=0}^{N-1} (x[n] - A -Bn)^{2} \bigg]
\end{aligned} }
\end{equation}
pdf에 기반한 (\ref{eq:intro1})과 같은 추정은 고전적인 추정 방법으로써 파라미터 $\theta$를 우리가 모르지만 고정된 상수로 보는 빈도주의(frequentist) 관점으로 해석될 수 있다. 이와 달리 현대적인 추정 방법은 파라미터 $\theta$ 또한 별도의 확률 변수로 해석하는 베이지안(bayesian) 관점을 주로 사용한다.
\begin{equation}
\begin{aligned}
& \text{Frequentist: } \quad \underbrace{x[n]}_{\text{r.v.}} = \underbrace{\theta}_{\text{deterministic}}+ w[n] \\
& \text{Bayesian: } \quad \underbrace{x[n]}_{\text{r.v.}} = \underbrace{\theta}_{\text{r.v.}} + w[n] \\
\end{aligned}
\end{equation}
베이지안 관점에서는 데이터 $\mathbf{x}$와 파라미터 $\theta$가 둘 다 확률변수이므로 둘 사이의 결합 확률 분포(joint pdf)을 사용하여 확를을 표현한다.
\begin{equation}
\begin{aligned}
p(\mathbf{x}, \theta) = p(\mathbf{x} | \theta) p(\theta)
\end{aligned}
\end{equation}
- $p(\mathbf{x}, \theta)$: joint pdf
- $p(\mathbf{x} | \theta)$ : conditional pdf
- $p(\theta)$: prior pdf
베이지안 관점에 대한 설명은 챕터 10에서 자세히 다룬다.
1.2. Assessing Estimator Performance
위와 같이 100일간 측정한 노이즈를 포함한 몸무게 데이터가 주어졌다고 가정하자. 위 데이터는 아래와 같이 모델링할 수 있다.
\begin{equation}
\begin{aligned}
x[n] = A + w[n]
\end{aligned}
\end{equation}
- $w[n]$ : 평균이 0인 노이즈
- $x[n]$ : 측정된 데이터
- $A$ : 추정하고자 하는 파라미터
일반적으로 $A$를 다음과 같이 데이터들의 평균으로 추정하는 것이 합리적일 것이다.
\begin{equation}
\boxed{ \begin{aligned}
\hat{A} = \frac{1}{N} \sum_{n=0}^{N-1} x[n]
\end{aligned} }
\end{equation}
- $\hat{A}$: A의 추정값 1
여기에 다음과 같은 질문을 할 수 있다.
- $\hat{A}$는 실제 $A$와 얼마나 가까울까?
- 평균말고 더 좋은 추정 방법은 없을까?
다음과 같이 다른 추정 방법을 사용하여 $A$를 추정할 수 있다.
\begin{equation}
\boxed{ \begin{aligned}
\breve{A} = x[0]
\end{aligned} }
\end{equation}
- $\breve{A}$: A의 추정값 2
직관적으로 우리는 모든 데이터(또는 무한개의 데이터)를 관측하는게 아닌 이상 좋은 추정값을 얻는 것이 어렵다는 것을 알 수 있다. 실제로 $\breve{A} = 69.8$이고 $\hat{A} = 71.1$이어서 $\breve{A}$가 $A=70$에 더 가깝다. 이런 경우 $\breve{A}$가 더 좋은 추정값이라고 볼 수 있을까? 당연히 아니다.
추정값(estimator) $\hat{A}$는 확률변수 $x[n]$에 대한 함수이므로 $\hat{A}$ 역시 확률변수가 된다. 따라서 추정값 역시 노이즈로 인해 다양한 결과물들을 도출할 수 있다. $\breve{A}$가 $A$에 더 가깝다는 사실은 주어진 $x[n]$의 예제에 대해서만을 의미한다. 따라서 추정값의 성능을 평가하기 위해서는 반드시 통계적으로 접근해야 한다. 예를 들어 여러번의 실험을 통해 데이터를 수집하고 이를 반복적으로 추정하는 방법이 존재한다.
데이터를 여러번 수집한 후 추정값 $\hat{A}, \breve{A}$의 기대값(expectation)을 계산하면 다음과 같다.
\begin{equation}
\begin{aligned}
\mathbb{E}(\hat{A}) & = \mathbb{E}\bigg ( \frac{1}{N} \sum_{n=0}^{N-1} x[n] \bigg) \\
& = \frac{1}{N} \sum_{n=0}^{N-1} \mathbb{E}( x[n] ) \\
& = \frac{1}{N} \sum_{n=0}^{N-1} A \\
& = A
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
\mathbb{E}(\breve{A}) & = \mathbb{E}(x[0]) \\
& = A
\end{aligned}
\end{equation}
따라서 두 추정값의 성능은 동일한 것일까? $\hat{A}$가 $\breve{A}$보다 더 좋은 추정값임을 증명하기 위해서는 추정의 분산이 더 작음을 입증해야 한다.
\begin{equation}
\begin{aligned}
\text{var}(\hat{A}) & = \text{var} \bigg( \frac{1}{N} \sum_{n=0}^{N-1} x[n] \bigg) \\
& = \frac{1}{N^{2}} \sum_{n=0}^{N-1} \text{var}(x[n]) \\
& =\frac{1}{N^{2}} N \sigma^{2} \\
& = \frac{\sigma^{2}}{N}
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
\text{var}(\breve{A}) & = \text{var}(x[0]) \\
& = \sigma^{2} \\
& > \text{var}(\hat{A})
\end{aligned}
\end{equation}
따라서 이를 통해 얻을 수 있는 결론은 다음과 같다.
- 추정값(estimator)는 확률변수(random variable)이다. 따라서 추정값의 성능은 반드시 통계적 방법이나 pdf를 통해 판단되어야 한다.
- 컴퓨터 시뮬레이션을 사용하여 추정값을 평가하는 방법은 파라미터에 대한 통찰을 얻기엔 충분히 좋지만 이를 절대적인 값으로 해석하면 안된다. 운이 좋은 경우 추정값은 소수점 오차를 가진 정확도로 구할 수 있지만 운이 나쁜 경우에는(데이터가 부족하거나 에러 값이 들어 있거나) 잘못된 추정값을 얻을 수 있다.
2. Minimum Variance Unbiased Estimation
본 섹션에서 나오는 추정값들은 고전적인 빈도주의 관점에서 파라미터 $\theta$가 고정된 값으로 주어졌다고 가정한다.
2.1. Unbiased Estimators
추정값이 불편성(unbiased)을 지닌다는 의미는 추정값의 평균이 파라미터의 참 값과 동일하다는 의미와 동치이다. 일반적으로 파리미터는 특정 범위 $a < \theta < b$ 안에 존재하므로 불편추정값(unbiased estimator, 또는 불편추정량)이란 다음과 같이 수학적으로 정의할 수 있다.
\begin{equation}
\begin{aligned}
\mathbb{E}(\hat{\theta}) = \theta
\end{aligned}
\end{equation}
만약 추정값이 편향(biased)되어 있다면 $\mathbb{E}(\hat{\theta}) \neq \theta$이고 편향 $b(\theta)$은 다음과 같이 계산할 수 있다.
\begin{equation}
\boxed{ \begin{aligned}
b(\theta) = \mathbb{E}(\hat{\theta}) - \theta
\end{aligned} }
\end{equation}
2.1.1. Example 2.1 and Example 2.2
다음과 같은 데이터가 주어졌다고 하자.
\begin{equation}
\begin{aligned}
x[n] = A + w[n] \quad n =0,1\cdots,N-1
\end{aligned}
\end{equation}
- $A$: 추정하고자 하는 파라미터, $-\infty < A < \infty$
- $w[n]$: WGN
이에 대한 일반적인 추정값 $\hat{A}$는 다음과 같이 데이터의 평균으로 예측할 수 있다.
\begin{equation}
\begin{aligned}
\hat{A} = \frac{1}{N} \sum_{n=0}^{N-1} x[n]
\end{aligned}
\end{equation}
선형성에 의해 기대값(expectation)은 다음과 같이 정의된다.
\begin{equation}
\begin{aligned}
\mathbb{E}(\hat{A}) & = \mathbb{E}\bigg ( \frac{1}{N} \sum_{n=0}^{N-1} x[n] \bigg) \\
& = \frac{1}{N} \sum_{n=0}^{N-1} \mathbb{E}( x[n] ) \\
& = \frac{1}{N} \sum_{n=0}^{N-1} A \\
& = A
\end{aligned}
\end{equation}
따라서 평균 추정값 $\hat{A}$는 불편추정값(unbiased estimator)이다. 만약 다음과 같은 추정값 $\breve{A}$가 있다고 가정해보자.
\begin{equation}
\begin{aligned}
\breve{A} = \frac{1}{2N} \sum_{n=0}^{N-1} x[n]
\end{aligned}
\end{equation}
$\breve{A}$의 기대값은 다음과 같다.
\begin{equation}
\begin{aligned}
\mathbb{E}(\breve{A}) & = \frac{1}{2}A \\
& = A \quad \text{ if } A = 0 \\
& \neq A \quad \text{ if } A \neq 0
\end{aligned}
\end{equation}
따라서 $\breve{A}$는 편의추정값(biased estimator)이다.
불편추정값이 반드시 좋은 추정값을 의미할까? 추정값이 불편성을 지닌다고 해서 반드시 좋은 추정값이라는 의미는 아니다. 불편성의 의미는 오직 추정값의 기대값(expectation)이 실제 값과 동일하다는 의미만 지닌다. 이와 반대로 편의추정값은 시스템의 노이즈를 포함하여 모델링한 값일 수 있다. 하지만 영구적인 편향성은 항상 안 좋은 추정값을 가진다. 예를 들면 다음과 같이 동일 파라미터 $\theta$에 대한 여러 추정값 $\{ \hat{\theta}_{1}, \hat{\theta}_{2}, \cdots, \hat{\theta}_{n}\}$이 주어졌을 때 가장 합리적인 방법은 이들의 기대값을 구하는 것이다.
\begin{equation}
\begin{aligned}
\hat{\theta} = \frac{1}{n} \sum _{i=1}^{n} \hat{\theta}_{i}
\end{aligned}
\end{equation}
만약 모든 추정값들이 불편성을 지니고 서로 독립이라면 다음 공식이 성립한다.
\begin{equation}
\begin{aligned}
\mathbb{E}(\hat{\theta}) = \theta
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
\text{var}(\hat{\theta}) & = \frac{1}{n^2} \sum_{i=1}^{n} \text{var} (\hat{\theta}_{i}) \\
& = \frac{\text{var}(\hat{\theta}_{1})}{n}
\end{aligned}
\end{equation}
따라서 위 식에서 보다시피 많은 수($n \uparrow$)의 추정값을 사용할 수록 분산은 작아진다. 만약 $n \rightarrow \infty$이면 분산은 0이 되고 $\hat{\theta} \rightarrow \theta$가 된다. 하지만 편의추정량의 경우 $\mathbb{E}(\hat{\theta}_{i}) = \theta + b(\theta)$이므로 다음과 같은 기대값을 가진다.
\begin{equation}
\begin{aligned}
\mathbb{E}(\hat{\theta}) & = \frac{1}{n} \sum_{i=1}^{n} \mathbb{E} (\hat{\theta}_{i}) \\
& = \theta + b(\theta)
\end{aligned}
\end{equation}
따라서 $n$이 충분히 많다고 하더라도 편향 $b(\theta)$ 값은 제거되지 않으므로 실제 추정값으로 수렴하지 않는다.
2.2. Minimum Variance Criterion
최적의 추정값을 찾기 위해서는 최적의 criterion을 사용해야 한다. 널리 사용되는 criterion 중 하나가 mean square error(MSE)이다.
\begin{equation}
\boxed{ \begin{aligned}
\text{mse}(\hat{\theta}) = \mathbb{E}\bigg[ (\hat{\theta} - \theta)^{2} \bigg]
\end{aligned} }
\end{equation}
- $\theta$ : 추정해야 할 파라미터. 빈도주의 관점에 의해 $\theta$는 고정된 상수 값을 의미한다. 즉, 확률변수가 아니다.
MSE는 실제 값 $\theta$과 추정값 $\hat{\theta}$의 평균 제곱 편차를 측정한다. MSE는 널리 사용되는 criterion 중 하나이지만 아쉽게도 편향에 의한 영향을 받는다. 위 식에 $\pm \mathbb{E}(\hat{\theta})$를 추가한 후 식을 전개하면 다음과 같다.
\begin{equation}
\begin{aligned}
\text{mse}(\hat{\theta}) & = \mathbb{E}\bigg\{ \Big[ \Big( \hat{\theta}- \mathbb{E}(\hat{\theta}) \Big) + \Big( \mathbb{E}(\hat{\theta}) - \theta \Big) \Big]^{2} \bigg\} \\
& = \text{var}(\hat{\theta}) + \bigg[ \mathbb{E}(\hat{\theta}) - \theta \bigg]^2 \\
& = \text{var}(\hat{\theta}) + b^{2}(\theta)
\end{aligned}
\end{equation}
따라서 최적의 추정값을 찾기 위한 criterion을 고려할 때 MSE를 최소화해주는 minimum MSE(MMSE) 추정값을 고려하면 안된다. MMSE의 대안으로 편향이 0이고 분산을 최소화하는 추정값을 사용해야 하는데 이를 minimum variance unbiased estimator(MVUE)라고 한다. 불편추정값의 분산을 최소화하는 과정은 pdf 관점에서 $p(\hat{\theta} - \theta)$를 0 주변에 집중시키는 효과가 있다. 따라서 추정 오차가 커질 가능성이 작아진다.
2.3. Existence of the Minimum Variance Unbiased Estimator
독자는 MVUE가 실제로 존재하는지 여부에 대해 궁금증이 생길 수 있다. 즉, 모든 파라미터 $\theta$에 대해 최소의 분산을 가지며 불편향된 추정값이 존재하는지 궁금할 수 있다. 결론만 말하자면 MVUE는 항상 존재하는 것은 아니다.
2.4. Finding the Minimum Variance Unbiased Estimator
만약 MVUE가 존재한다고해도 이를 찾는 것이 불가능할 수 있다. MVUE를 찾을 수 있는 절대적인 방법이란 아직 알려지지 않았다. 하지만 이를 가능하게 해주는 몇몇 기준들은 존재한다.
- Cramer-Rao lower bound(CRLB)를 결정하고 추정값들이 이를 만족하는지 확인한다.
- Rao-Blackwell-Lehmann-Scheffe(RBLS) 이론을 적용한다.
- 추정값이 불편성(unbiased) 뿐만아니라 선형(linear) 특성이 있다는 제약조건 하에 분산을 최소화하는 MVUE를 구한다.
1,2 방법을 사용하면 MVUE를 구할 수 있고 3 방법은 MVUE가 데이터에 대하여 선형인 경우에만 적용된다.
- 1,2) CRLB는 임의의 불편추정값에 대하여 분산의 하한선(lower bound)를 결정하게 해준다. 만약 어떤 추정값이 CRLB와 동일한 분산 값을 가진다면 이는 반드시 MVUE가 된다. CRLB와 동일한 분산 값을 가지지 않는다고 하더라도 MVUE가 존재할 수 있다. 이럴 때는 RBLS를 적용한다. RBLS는 충분통계량(sufficient statistics)을 먼저 구한 후 충분통계량에 대한 추정을 수행하는데 이 때 추정값이 $\theta$에 대한 불편추정값이 된다.
- 3) 이는 추정값이 선형이어야 하는 제약조건을 가진다. 이는 때때로 강력한 제약조건이지만 최적의 선형 추정값을 구할 수 있다.
3. Cramer-Rao Lower Bound
어떠한 불편추정값(unbiased estimator)의 분산의 하한선(lower bound)를 결정할 수 있다는 것은 실제 추정 문제에서 매우 유용하게 사용된다. 최선의 경우 특정 추정값이 MVUE임을 바로 구할수도 있다. 그렇지 않은 경우라 하더라도 여러 불편추정량에 대한 벤치마크 용도로 활용될 수도 있다. 이는 신호 추정 분야에서 매우 유용하게 활용된다. Cramer-Rao lower bound(CRLB) 이외에도 [McAulay and Hofstetter 1971, Kendall and Stuart 1979, Seidman 1970, Ziv and Zakai 1969]와 같이 분산의 한계(bound)를 결정하는 알고리즘들이 존재하지만 CRLB가 이들 중 가장 쉽게 한계를 구할 수 있다.
3.1. Estimator Accuracy Considerations
추정에 사용되는 정보는 일반적으로 관측된 데이터로부터 얻을 수 있고 관측 데이터는 노이즈를 포함하기 때문에 일반적으로 pdf로 표현될 수 있다. 따라서 추정의 정확도는 당연히 pdf와 직접적인 관련이 있다. 만약 파라미터가 pdf에 영향을 거의 주지 않는 최악의 경우에는 좋은 추정값을 얻는 것이 어려울 것이다. 따라서 파라미터가 pdf에 영향을 많이 줄수록 추정의 정확도는 올라간다.
3.1.1. Example 3.1 - PDF Dependence on Unknown Parameter
만약 하나의 데이터가 샘플링되었다고 가정해보자
\begin{equation}
\begin{aligned}
x[0] = A + w[0]
\end{aligned}
\end{equation}
이 때, $w[n]$는 $\mathcal{N}(0, \sigma^{2})$을 따르는 white gaussian noise(WGN)이다. 일반적으로 좋은 추정값(estimator)이란 $\sigma^{2}$가 작은 추정값임을 알 수 있다. 그리고 좋은 불편추정값은 $\hat{A} = x[0]$임을 알 수 있다.
분산 $\sigma^{2}$ 값이 작을 수록 좋은 추정값임을 설명하기 위해 아래와 같은 예제를 들 수 있다. 만약 서로 다른 두 분산 $\sigma_{1}=1/3$와 $\sigma_{2}=1$이 주어졌고 $x[0]=3$인 경우를 가정해보자. 이에 대한 pdf는 다음과 같다.
\begin{equation}
\begin{aligned}
p_{i}(x[0];A) = \frac{1}{ \sqrt{ 2\pi\sigma_{i}^{2} }} \exp\bigg[ -\frac{1}{2\sigma_{i}^{2}} (x[0] - A)^{2} \bigg] \quad i=1,2
\end{aligned}
\end{equation}
$\sigma_{1}^2 < \sigma_{2}^2$이므로 우리는 $p_{1}(\cdot)$이 $A$를 더 잘 추정하고 있다고 판단할 수 있다. 예를 들어 $A=4.5$일 확률은 $p_{1}$보다 $p_{2}$가 더 높으므로 우리는 $p_{1}$이 더 좁은 범위에 대한 확실한 확률분포를 가짐을 알 수 있다.
Likelihood Function v.s. Probability Density Function
만약 pdf $p(x; \theta)$가 $x$는 고정된 값이면서 동시에 파라미터 $\theta$에 대한 함수라면 이를 일반적으로 가능도함수(likelihood function)라고 부른다. 이와 반대로, $p(x;\theta)$가 $\theta$는 고정된 값이면서 동시에 $x$에 대한 함수라면 이는 일반적인 확률밀도함수(probability density function, pdf)라고 부른다. 두 개념을 비교한 그림은 아래와 같다. $p(x;\theta)$는 다음과 같다.
\begin{equation} \begin{aligned} p(x;\theta) = \frac{1}{ \sqrt{ 2\pi\sigma^{2} }} \exp\bigg[ -\frac{1}{2\sigma^{2}} (x - \theta)^{2} \bigg] \end{aligned} \end{equation}
그림에서 보는 것과 같이 pdf와 가능도함수는 모양만 같을 뿐 이를 해석하는 관점이 서로 다르다. pdf의 경우 파라미터가 주어졌을 때 특정 구간 내에서 확률을 구하는 것에 관심이 있다면 가능도함수는 데이터가 주어졌을 때 이를 가장 잘 설명하는 파라미터는 무엇인가? 에 관심이 있다.
지금까지 설명한 pdf는 전부 데이터 $x$가 주어졌을 때 파라미터 $A$를 찾는 형태이기 때문에 이는 가능도함수(likelihood function) 관점에서 해석할 수 있다. 가능도함수의 뾰족한 정도(sharpness)는 추정값이 얼마나 정확한 지에 대한 정확도를 판단하는데 사용된다. 수학적으로 곡선의 뾰족한 정도는 함수의 2차 미분을 수행하여 곡률(curvature)를 구함으로써 얻을 수 있다. 이 때, pdf 특성 상 exponential 항이 존재하여 계산이 어렵기 때문에 일반적으로 log를 취한 값을 사용하는데 이를 로그 가능도함수(log likelihood function)이라고 한다.
\begin{equation}
\begin{aligned}
\ln p(x[0];A) = -\frac{1}{2} \ln(2\pi\sigma^{2}) -\frac{1}{2 \sigma^2} (x[0] - A)^{2}
\end{aligned}
\end{equation}
파라미터 $A$에 대한 1차 미분을 수행하면 다음과 같다.
\begin{equation}
\begin{aligned}
\frac{\partial \ln p(x[0];A)}{\partial A} = \frac{1}{\sigma^2} (x[0] - A)
\end{aligned}
\end{equation}
다시 한번 2차 미분을 취한 후 양변에 음수를 곱하면 다음과 같다.
\begin{equation} \label{eq:crlb1}
\begin{aligned}
- \frac{\partial^{2} \ln p(x[0];A)}{\partial A^{2}} = \frac{1}{\sigma^{2}}
\end{aligned}
\end{equation}
위 식의 의미는 $\sigma^{2}$가 감소할수록 원래 함수의 곡률(curvature)을 의미하는 $-\frac{\partial^{2} \ln p(x[0];A)}{\partial A^{2}}$는 증가하는 것을 의미한다. 앞서 말한 추정값 $\hat{A} = x[0]$의 분산은 다음과 같다.
\begin{equation}
\begin{aligned}
\text{var}(\hat{A}) = \sigma^{2}
\end{aligned}
\end{equation}
여기에 (\ref{eq:crlb1})을 대입하면 다음과 같다.
\begin{equation}
\begin{aligned}
\text{var}(\hat{A}) = \frac{1}{- \frac{\partial^{2} \ln p(x[0];A)}{\partial A^{2}}}
\end{aligned}
\end{equation}
위와 같은 간단한 예제에서는 로그 가능도함수의 2차 미분값이 데이터 $x[0]$에 독립이지만 일반적으로는 2차 미분값은 데이터 $\mathbf{x}=[x[0], \cdots, x[n]]$에 종속적이다. 곡률의 정확한 표현은 다음과 같이 나타낼 수 있다.
\begin{equation}
\boxed{ \begin{aligned}
-\mathbb{E} \bigg[ \frac{\partial^{2} \ln p(x[0];A)}{\partial A^{2}} \bigg]
\end{aligned} }
\end{equation}
위 식을 통해 다양한 관측 데이터가 주어졌을 때 로그 가능도함수의 평균적인 곡률을 측정할 수 있다.
3.2. Cramer-Rao Lower Bound
3.2.1. Theorem 3.1 (Cramer-Rao Lower Bound - Scalar Parameter)
pdf $p(\mathbf{x};\theta)$가 다음과 같은 정규 조건(regularity condition)을 만족하면
\begin{equation} \label{eq:crlb3}
\begin{aligned}
\mathbb{E} \bigg[ \frac{\partial \ln p(\mathbf{x};\theta)}{\partial \theta} \bigg] = 0 \quad \text{ for all } \theta
\end{aligned}
\end{equation}
임의의 불편추정값(unbiased estimator) $\hat{\theta}$의 분산은 다음 조건을 반드시 만족한다.
\begin{equation} \label{eq:crlb2}
\begin{aligned}
\text{var}(\hat{\theta}) \geq \frac{1}{ -\mathbb{E} \bigg[ \frac{\partial^{2} \ln p(\mathbf{x};\theta)}{\partial \theta^{2}} \bigg] }
\end{aligned}
\end{equation}
이 때 편미분은 실제 파라미터의 참 값 $\theta$에 대하여 수행되었다. 추가적으로 불편추정값의 분산이 하한선(lower bound)에 도달하려면 다음 조건을 반드시 만족해야 한다. (필요충분조건)
\begin{equation}
\begin{aligned}
\frac{\partial \ln p(\mathbf{x};\theta)}{\partial \theta} = I(\theta)(g(\mathbf{x}) - \theta)
\end{aligned}
\end{equation}
위 식은 로그 가능도함수의 1차 미분이 위와 같은 임의의 함수 $I, g$의 곱셈 형태로 표현되어야 한다는 뜻이다. 이 때 불편추정값은 $\hat{\theta} = g(\mathbf{x})$가 되며 $\hat{\theta}$는 MVUE를 만족한다. 이 때의 최소 분산 값은 $1/I(\theta)$이 되며 이를 Fisher information이라고 한다. Fisher information은 일반적으로 다음과 같이 정의한다.
\begin{equation} \label{eq:crlb-fisher}
\begin{aligned}
I(\theta) = -\mathbb{E} \bigg[ \frac{\partial^{2} \ln p(\mathbf{x};\theta)}{\partial \theta^{2}} \bigg]
\end{aligned}
\end{equation}
따라서 MVUE의 분산은 다음과 같이 표현하기도 한다.
\begin{equation} \label{eq:crlb-fisher2}
\begin{aligned}
\text{var}(\hat{\theta}) = \frac{1}{ I(\theta) } \quad \cdots \text{ for MVUE}
\end{aligned}
\end{equation}
(\ref{eq:crlb2})에서 2차 미분값은 $\mathbf{x}$에 종속적이기 때문에 기대값은 정의에 따라 다음과 같이 전개된다.
\begin{equation}
\begin{aligned}
\mathbb{E} \bigg[ \frac{\partial^{2} \ln p(\mathbf{x};\theta)}{\partial \theta^{2}} \bigg] = \int \frac{\partial^{2} \ln p(\mathbf{x};\theta)}{\partial \theta^{2}} p(\mathbf{x} ; \theta) d\mathbf{x}
\end{aligned}
\end{equation}
3.2.2. Example 3.3 - DC Level in White Gaussian Noise
다음과 같은 여러 관측 데이터들이 주어졌다고 하자.
\begin{equation}
\begin{aligned}
x[n] = A + w[n] \quad n=0,1,\cdots,N-1
\end{aligned}
\end{equation}
- $w[n] \sim \mathcal{N}(0, \sigma^{2})$ : WGN
파라미터 $A$에 대한 CRLB를 유도해보면 다음과 같다. 우선 모든 관측값 $\mathbf{x}$에 대한 가능도함수를 구한다.
\begin{equation}
\begin{aligned}
p(\mathbf{x};A) & = \Pi_{n=0}^{N-1} \frac{1}{ \sqrt{ 2\pi\sigma^{2} }} \exp\bigg[ -\frac{1}{2\sigma^{2}} (x[n] - A)^{2} \bigg] \\
& = \frac{1}{(2\pi\sigma^{2})^{\frac{N}{2}}} \exp\bigg[ -\frac{1}{2\sigma^{2}} \sum_{n=0}^{N-1}(x[n] - A)^{2} \bigg]
\end{aligned}
\end{equation}
로그 가능도 함수는 다음과 같다.
\begin{equation}
\begin{aligned}
\ln p(\mathbf{x};A) = -\ln[ (2\pi\sigma^{2})^{\frac{N}{2}} ] -\frac{1}{2\sigma^{2}} \sum_{n=0}^{N-1}(x[n] - A)^{2} \\
\end{aligned}
\end{equation}
1차 미분을 수행하면 다음과 같다.
\begin{equation}
\begin{aligned}
\frac{\partial \ln p(\mathbf{x};A)}{\partial A} & = \frac{\partial}{\partial A} \bigg[ -\ln[ (2\pi\sigma^{2})^{\frac{N}{2}} ] -\frac{1}{2\sigma^{2}} \sum_{n=0}^{N-1}(x[n] - A)^{2} \bigg] \\
& = \frac{1}{\sigma^{2}} \sum_{n=0}^{N-1}(x[n]-A) \\
& = \frac{N}{\sigma^{2}} (\bar{x} - A)
\end{aligned}
\end{equation}
- $\bar{x}$ : $\mathbf{x}$의 평균
2차 미분을 수행하면 다음과 같다.
\begin{equation}
\begin{aligned}
\frac{\partial^{2} \ln p(\mathbf{x};A)}{\partial A^{2}} & = -\frac{N}{\sigma^{2}}
\end{aligned}
\end{equation}
위 식은 파라미터가 없는 상수임에 유의한다. (\ref{eq:crlb2}) 식에 이를 대입하면 다음과 같다.
\begin{equation}
\begin{aligned}
\text{var}(\hat{A}) \geq \frac{\sigma^{2}}{N}
\end{aligned}
\end{equation}
위 식이 정규 조건을 만족하는가? (\ref{eq:crlb3})를 보면 적용해보면 만족하는 것을 알 수 있다.
\begin{equation}
\begin{aligned}
\mathbb{E} \bigg[ \frac{\partial \ln p(\mathbf{x};\theta)}{\partial \theta} \bigg] & =\mathbb{E} \bigg[ \frac{N}{\sigma^{2}} (\bar{x} - A) \bigg] = 0
\end{aligned}
\end{equation}
- $\mathbb{E}(\bar{x}) = A$
따라서 CRLB 정의에 따라 임의의 불편추정값 $\breve{A}$의 분산이 $\text{var}(\breve{A}) = \frac{\sigma^{2}}{N}$을 만족하면 이는 반드시 MVUE가 된다. 위 예제에서는 $\breve{A}=\bar{x}$가 MVUE가 된다. 그리고 이러한 추정값을 efficient하다라고 한다. 관측 데이터를 효율적(efficient)으로 사용하여 추정하였다는 의미이다.
3.3. Transformation of Parameters
실제 추정 문제에서는 우리가 추정하고자 하는 파라미터가 단순한 $\theta$가 아닌 $\theta$의 함수 형태를 추정해야 하는 일이 자주 발생한다. 이전 예제에서도 단순히 $A$를 추정하는 것이 아닌 $A^{2}$를 추정하고 싶을 수도 있다. 만약 A의 CRLB를 알고 있는 경우에는 $A^{2}$의 CRLB도 쉽게 구할 수 있고 $A$와 관련된 어떤 함수라도 구할 수 있다.
추정하고자 하는 파라미터가 $\alpha = g(\theta)$와 같이 $\theta$에 대한 함수일 경우 CRLB는 다음과 같다.
\begin{equation}
\begin{aligned}
\text{var}(\hat{\alpha}) \geq \frac{ \bigg(\frac{\partial g}{\partial \theta}\bigg)^{2} }{ -\mathbb{E} \bigg[ \frac{\partial^{2} \ln p(\mathbf{x};\theta)}{\partial \theta^{2}} \bigg] }
\end{aligned}
\end{equation}
만약 $\alpha = g(A) = A^{2}$인 경우 CRLB는 다음과 같다.
\begin{equation}
\begin{aligned}
\text{var}(\hat{\alpha}) \geq \frac{(2A)^2}{N/\sigma^{2}} = \frac{4A^{2}\sigma^{2}}{N}
\end{aligned}
\end{equation}
이전 Example 3.3 예제에서 $\hat{A} = \bar{x}$는 $A$에 대하여 efficient함을 보였다. 따라서 $\bar{x}^{2}$ 역시 $A^{2}$에 대하여 efiicient하다고 예상할 수 있으나 이는 사실이 아니다. 그 이전에 $\hat{A^{2}} = \bar{x}^{2}$는 $A^{2}$에 대한 불편추정값 조차 아니다. 즉, 편향(bias)이 존재한다.
\begin{equation}
\begin{aligned}
\mathbb{E}(\bar{x}^{2}) & = \mathbb{E}^{2}(\bar{x}) + \text{var}(\bar{x}) \\
& = A^{2} + \frac{\sigma^2}{N} \\
& \neq A^{2}
\end{aligned}
\end{equation}
따라서 우리는 위 예제를 통해 추정값의 efficiency는 비선형 변환에 의해 보존되지 않는 것을 알 수 있다. 하지만 linear 또는 affine 변환에 대하여는 efficiency가 보존됨을 쉽게 보일 수 있다.
만약 $\hat{\theta}$가 $\theta$에 대하여 efficient하고 $\alpha = g(\theta) = a\theta + b$와 같은 affine 변환이 주어진 경우를 확인해보자.
\begin{equation}
\begin{aligned}
\hat{\alpha} = g(\hat{\theta}) = a \hat{\theta}+b
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
\mathbb{E}(a \hat{\theta}+b) & = a\mathbb{E}(\hat{\theta}) + b \\
& = a \theta + b \\
& = g(\theta)
\end{aligned}
\end{equation}
\begin{equation} \label{eq:crlb4}
\begin{aligned}
\text{var}(\hat{\alpha}) & = \text{var}(a \hat{\theta} + b) \\
& = a^{2} \text{var}(\hat{\theta})
\end{aligned}
\end{equation}
$g(\theta)$의 CRLB를 보면 다음과 같다.
\begin{equation} \label{eq:crlb5}
\begin{aligned}
\text{var}(\hat{\alpha}) & \geq \frac{ \bigg(\frac{\partial g}{\partial \theta}\bigg)^{2} }{ I(\theta) } \\
& = \bigg(\frac{\partial g}{\partial \theta}\bigg)^{2} \text{var}(\hat{\theta}) \\
& = a^{2} \text{var}(\hat{\theta})
\end{aligned}
\end{equation}
위 식에서 (\ref{eq:crlb4}), (\ref{eq:crlb5})의 분산이 동일하기 때문에 $\hat{\alpha}$ 역시 MVUE이면서 동시에 efficient함을 알 수 있다.
앞서 보았듯이 efficiency는 linear 또는 affine 변환에서만 보존되는 것을 확인하였다. 하지만 데이터가 충분히 큰 경우에는 비선형 변환도 근사적으로(approximately) efficiency가 보존된다. 다시 이전 예제 $\alpha = g(A) = A^{2}$로 돌아가서 데이터 $x[n]$의 평균을 $\bar{x} = \frac{1}{N} \sum_{n=0}^{N-1} x[n]$이라고 할 때 $\bar{x}^{2}$는 $N$이 충분히 클 경우 근사적으로 편향이 제거된다.
\begin{equation} \label{eq:crlb6}
\begin{aligned}
\text{var}(\bar{x}^{2}) & = \mathbb{E}(\bar{x}^{4}) - \mathbb{E}^{2}(\bar{x}^{2}) \\
& = \frac{4A^{2}\sigma^{2}}{N} + \frac{2 \sigma^{2}}{N^2}\\
& \approx \frac{4A^{2}\sigma^{2}}{N} \quad \cdots \text{ for } N \rightarrow \infty
\end{aligned}
\end{equation}
$g(A) = A^{2}$의 CRLB는 다음과 같다.
\begin{equation} \label{eq:crlb7}
\begin{aligned}
\text{var}(\hat{\alpha}) & \geq \frac{(2A)^2}{N/\sigma^{2}} \\
& = \frac{4A^{2}\sigma^{2}}{N}
\end{aligned}
\end{equation}
(\ref{eq:crlb6}), (\ref{eq:crlb6})가 $N$이 충분히 큰 경우에 대하여 서로 동일하기 때문에 데이터가 많은 경우에는 비선형 변환에 대한 추정값도 MVUE가 되며 동시에 efficient함을 알 수 있다.
다른 방법을 사용하여 비선형 변환이 근사적으로 efficient함을 보일 수 있다. 확률분포 관점에서 봤을 때 $N$이 커질수록 $\hat{\mathbf{x}}$는 $A$ 주변으로 뾰족해지는 경향이 있다. 이에 따라 $\pm 3\sigma$ 사이의 간격은 좁아지게 되고 좁은 영역에 대해여 비선형 변환을 수행하면 근사적으로 선형 변환을 한 것과 유사한 효과를 얻는다. 이를 $\bar{x} = A$ 지점에서 테일러 1차 근사를 통해 표현하면 다음과 같다.
\begin{equation}
\begin{aligned}
g(\bar{x}) \approx g(A) + \frac{dg(A)}{dA} (\bar{x} - A)
\end{aligned}
\end{equation}
이러한 경우를 점근적으로(asymptotically) efficient하다고 한다. 이 때 기대값은 다음과 같다.
\begin{equation}
\begin{aligned}
\mathbb{E}[g(\bar{x})] = g(A) = A^2
\end{aligned}
\end{equation}
분산은 다음과 같다.
\begin{equation}
\begin{aligned}
\text{var}[g(\bar{x})] & = \bigg[ \frac{dg(A)}{dA} \bigg]^{2} \text{var}(\bar{x}) \\
& = \frac{(2A)^2 \sigma^2 }{N} \\
& = \frac{4A^2 \sigma^2}{N}
\end{aligned}
\end{equation}
즉 추정값은 CRLB에 점근적으로(asymptotically) 도달하는 것을 알 수 있다. 따라서 비선형 변환은 점근적으로 efficient하다.
3.4. Extension to a Vector Parameter
지금까지는 추정하려는 파라미터 $\theta$가 스칼라 값이었다. 해당 섹션에서는 파라미터가 여러개인 벡터 파리미터 $\boldsymbol{\theta} = [\theta_{1}, \theta_{2},\cdots, \theta_{p}]^{\intercal}$ 케이스에 대해 다룬다. $\hat{\boldsymbol{\theta}}$는 $\boldsymbol{\theta}$에 대한 불편추정값(unbiased estimator)이라고 가정한다. 벡터 파라미터의 각 원소의 분산은 다음과 같이 나타낼 수 있다.
\begin{equation} \label{eq:crlb8}
\begin{aligned}
\text{var}(\hat{\theta}_{i}) \geq [\mathbf{I}^{-1}(\boldsymbol{\theta})]_{ii}
\end{aligned}
\end{equation}
- $\mathbf{I}(\boldsymbol{\theta}) \in \mathbb{R}^{p \times p}$ : Fisher information 행렬
이는 (\ref{eq:crlb-fisher2})의 벡터 버전임을 알 수 있다. 일반적으로 스칼라 버전에서 행렬은 벡터 버전에서 역행렬로 표현된다. Fisher information 행렬을 자세히 나타내면 다음과 같다.
\begin{equation}
\begin{aligned}
\big[ \mathbf{I}(\boldsymbol{\theta}) \big]_{ij} = -\mathbb{E} \bigg[ \frac{\partial^{2} \ln p(\mathbf{x};\theta)}{\partial \theta_{i} \partial \theta_{j}} \bigg]
\end{aligned}
\end{equation}
- $i = 1,2,\cdots, p$
- $j = 1,2,\cdots, p$
$p=1$인 스칼라 케이스는 $\mathbf{I}(\boldsymbol{\theta}) = I(\theta)$가 된다. 스칼라 버전과 동일하게 불편추정값의 분산이 하한선(lower bound)에 도달하려면 다음 조건을 반드시 만족해야 한다. (필요충분조건)
\begin{equation}
\begin{aligned}
\frac{\partial \ln p(\mathbf{x};\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = \mathbf{I}(\boldsymbol{\theta})(\mathbf{g}(\mathbf{x}) - \boldsymbol{\theta})
\end{aligned}
\end{equation}
위 식은 로그 가능도함수의 1차 미분이 위와 같은 임의의 함수 $\mathbf{I}, \mathbf{g}$의 곱셈 형태로 표현되어야 한다는 뜻이다. 이 때 불편추정값은 $\hat{\boldsymbol{\theta}} = \mathbf{g}(\mathbf{x})$가 되며 $\hat{\boldsymbol{\theta}}$는 MVUE를 만족한다. 이 때의 최소 분산 값은 $1/\mathbf{I}(\boldsymbol{\theta})$이 된다.
\begin{equation}
\begin{aligned}
\text{var}(\hat{ \boldsymbol{\theta}}) = 1/\mathbf{I}(\boldsymbol{\theta}) \quad \cdots \text{ for MVUE}
\end{aligned}
\end{equation}
3.4.1. Example 3.6 - DC Level in White Gaussian Noise (Revisited)
예제 3.3과 같이 관측 데이터들이 주어졌다고 하자.
\begin{equation}
\begin{aligned}
x[n] = A + w[n] \quad n=0,1,\cdots,N-1
\end{aligned}
\end{equation}
- $w[n] \sim \mathcal{N}(0, \sigma^{2})$ : WGN
이 때, 추정하고자 하는 파라미터가 벡터 파라미터 $ \boldsymbol{\theta} = [A, \sigma^{2}]^{\intercal}$인 경우를 생각해보자. 즉 $p=2$이다. Fisher information 행렬은 다음과 같이 나타낼 수 있다.
\begin{equation}
\begin{aligned}
\mathbf{I}(\boldsymbol{\theta}) = \begin{bmatrix}
-\mathbb{E} \bigg[ \frac{\partial^{2} \ln p(\mathbf{x};\theta)}{\partial A^2} \bigg] & -\mathbb{E} \bigg[ \frac{\partial^{2} \ln p(\mathbf{x};\theta)}{\partial A \partial \sigma^2} \bigg] \\
-\mathbb{E} \bigg[ \frac{\partial^{2} \ln p(\mathbf{x};\theta)}{\partial \sigma^2 \partial A} \bigg] &
-\mathbb{E} \bigg[ \frac{\partial^{2} \ln p(\mathbf{x};\theta)}{\partial {\sigma^{2}}^{2} } \bigg]
\end{bmatrix} \in \mathbb{R}^{2\times2}
\end{aligned}
\end{equation}
Fisher information 행렬은 대칭이며 동시에 positive definite한 특징을 가지고 있다. 예제 3.3의 로그 가능도함수는 다음과 같다.
\begin{equation}
\begin{aligned}
\ln p(\mathbf{x};\boldsymbol{\theta}) = -\frac{N}{2} \ln 2\pi - \frac{N}{2} \ln \sigma^{2} - \frac{1}{2\sigma^{2}} \sum_{n=0}^{N-1}(x[n] - A)^{2} \\
\end{aligned}
\end{equation}
로그 가능도함수의 1,2차 편미분은 다음과 같이 쉽게 구할 수 있다.
\begin{equation}
\begin{aligned}
\frac{ \ln p(\mathbf{x};\boldsymbol{\theta}) }{\partial A} & = \frac{1}{\sigma^{2}} \sum_{n=0}^{N-1}(x[n] - A) \\
\frac{ \ln p(\mathbf{x};\boldsymbol{\theta}) }{\partial \sigma^2} & = - \frac{N}{2 \sigma^2 } + \frac{1}{2\sigma^{4}} \sum_{n=0}^{N-1}(x[n] - A)^{2} \\
\frac{ \ln^2 p(\mathbf{x};\boldsymbol{\theta}) }{\partial A^2} & = -\frac{N}{\sigma^{2}} \\
\frac{ \ln^2 p(\mathbf{x};\boldsymbol{\theta}) }{\partial A \partial \sigma^{2}} & = -\frac{1}{\sigma^{4}} \sum_{n=0}^{N-1}(x[n] - A) \\
\frac{ \ln^2 p(\mathbf{x};\boldsymbol{\theta}) }{\partial {\sigma^{2}}^2} & = \frac{N}{2\sigma^4} - \frac{1}{\sigma^6} \sum_{n=0}^{N-1}(x[n] - A)^{2}
\end{aligned}
\end{equation}
이를 Fisher information 행렬에 대입하면 다음과 같다.
\begin{equation}
\begin{aligned}
\mathbf{I}(\boldsymbol{\theta}) = \begin{bmatrix}
\frac{N}{\sigma^2} & 0 \\ 0 & \frac{N}{2\sigma^4}
\end{bmatrix}
\end{aligned}
\end{equation}
흔한 경우는 아니지만 예제의 케이스는 역행렬을 매우 쉽게 구할 수 있다. 단순히 역수를 취해줌으로써 역행렬을 구하면 (\ref{eq:crlb8})와 같이 CRLB를 구할 수 있다.
\begin{equation}
\begin{aligned}
\text{var}(\hat{A})& \geq \frac{\sigma^2}{N} \\
\text{var}(\hat{\sigma^{2}}) & \geq \frac{2\sigma^4}{N}
\end{aligned}
\end{equation}
이 중 $\text{var}(\hat{A})$는 스칼라 케이스에서 $\sigma^{2}$ 값이 이미 주어진 경우와 동일한 CRLB를 가지는 것을 알 수 있다. 다시 말하자면 이러한 예제는 일반적인 상황에서는 참이 아니지만 예제의 경우 참이다.
3.4.2. Example 3.7 - Line Fitting
다음과 같은 line fitting 문제가 주어졌다고 가정해보자
\begin{equation}
\begin{aligned}
x[n] = A + Bn + w[n] \quad n=0,1,\cdots,N-1
\end{aligned}
\end{equation}
- $w[n] \sim \mathcal{N}(0, \sigma^{2})$ : WGN
이 때, y절편의 값 $A$와 기울기 $B$의 값을 찾아야 한다. 추정하고자 하는 파라미터는 $\theta=[A, \ B]^{\intercal}$이다. $p=2$ 케이스이므로 Fisher information 행렬은 다음과 같이 나타낼 수 있다.
\begin{equation}
\begin{aligned}
\mathbf{I}(\boldsymbol{\theta}) = \begin{bmatrix}
-\mathbb{E} \bigg[ \frac{\partial^{2} \ln p(\mathbf{x};\theta)}{\partial A^2} \bigg] & -\mathbb{E} \bigg[ \frac{\partial^{2} \ln p(\mathbf{x};\theta)}{\partial A \partial B^2} \bigg] \\
-\mathbb{E} \bigg[ \frac{\partial^{2} \ln p(\mathbf{x};\theta)}{\partial B \partial A} \bigg] &
-\mathbb{E} \bigg[ \frac{\partial^{2} \ln p(\mathbf{x};\theta)}{\partial B^{2} } \bigg]
\end{bmatrix} \in \mathbb{R}^{2\times2}
\end{aligned}
\end{equation}
가능도함수는 다음과 같다.
\begin{equation}
\begin{aligned}
p(\mathbf{x},\boldsymbol{\theta}) = \frac{1}{(2\pi\sigma^{2})^{\frac{N}{2}}} \exp\bigg[ -\frac{1}{2\sigma^{2}} \sum_{n=0}^{N-1}(x[n] - A- Bn)^{2} \bigg]
\end{aligned}
\end{equation}
로그 가능도함수의 1,2차 미분은 다음과 같다.
\begin{equation}
\begin{aligned}
\frac{ \ln p(\mathbf{x};\boldsymbol{\theta}) }{\partial A} & = \frac{1}{\sigma^{2}} \sum_{n=0}^{N-1}(x[n] - A - Bn) \\
\frac{ \ln p(\mathbf{x};\boldsymbol{\theta}) }{\partial B} & = \frac{1}{\sigma^{2}} \sum_{n=0}^{N-1}(x[n] - A - Bn)n \\
\frac{ \ln^2 p(\mathbf{x};\boldsymbol{\theta}) }{\partial A^2} & = -\frac{N}{\sigma^{2}} \\
\frac{ \ln^2 p(\mathbf{x};\boldsymbol{\theta}) }{\partial A \partial B} & = -\frac{1}{\sigma^{2}} \sum_{n=0}^{N-1}n \\
\frac{ \ln^2 p(\mathbf{x};\boldsymbol{\theta}) }{\partial {\sigma^{2}}^2} & = -\frac{1}{\sigma^{2}} \sum_{n=0}^{N-1}n^2
\end{aligned}
\end{equation}
이를 Fisher information 행렬에 대입하면 다음과 같다.
\begin{equation}
\begin{aligned}
\mathbf{I}(\boldsymbol{\theta}) & = \frac{1}{\sigma^2} \begin{bmatrix} N & \sum_{n=0}^{N-1}n \\
\sum_{n=0}^{N-1}n & \sum_{n=0}^{N-1}n^2
\end{bmatrix} \\
&= \frac{1}{\sigma^2} \begin{bmatrix}
N & \frac{N(N-1)}{2} \\
\frac{N(N-1)}{2} & \frac{N(N-1)(2N-1)}{6}
\end{bmatrix} \\
\end{aligned}
\end{equation}
$\mathbf{I}(\boldsymbol{\theta})$의 역행렬을 구해보면 다음과 같다.
\begin{equation}
\begin{aligned}
\mathbf{I}(\boldsymbol{\theta})^{-1} &=
\sigma^2 \begin{bmatrix}
\frac{2(2N-1)}{N(N+1)} & -\frac{6}{N(N+1)} \\
-\frac{6}{N(N+1)} & \frac{12}{N(N^2-1)}
\end{bmatrix} \\
\end{aligned}
\end{equation}
(\ref{eq:crlb8})와 같이 CRLB를 구해보면 다음과 같다.
\begin{equation}
\begin{aligned}
\text{var}(\hat{A})& \geq \frac{2(2N-1)\sigma^2}{N(N+1)} \\
\text{var}(\hat{B}) & \geq \frac{12\sigma^2}{N(N^2-1)}
\end{aligned}
\end{equation}
위 예제에서 보다시피 스칼라 파라미터 $\theta=A$만 추정했을 때와는 달리 $\boldsymbol{\theta}= [A, \ B]^{\intercal}$처럼 벡터 파라미터를 추정하면 CRLB는 커지는 것을 알 수 있다. 스칼라 파라미터만 추정했을 때 CRLB는 다음과 같다.
\begin{equation}
\begin{aligned}
\text{var}(\hat{A})& \geq \frac{\sigma^2}{N} \\
\end{aligned}
\end{equation}
따라서 $N \geq 2$인 경우에는 벡터 파라미터의 CRLB가 항상 스칼라 파라미터의 CRLB보다 크다.
\begin{equation}
\begin{aligned}
\frac{2(2N-1)\sigma^2}{N(N+1)} > \frac{\sigma^2}{N} \quad \cdots \text{ for } N \geq 2
\end{aligned}
\end{equation}
또한 특정 파라미터는 다른 파라미터보다 데이터 개수 $N$에 민감하게 반응할 수 있다.
\begin{equation}
\begin{aligned}
\frac{\text{CRLB}(\hat{A})}{\text{CRLB}(\hat{B})} = \frac{(2N-1)(N-1)}{6} > 1 \quad \cdots \text{ for } N \geq 3
\end{aligned}
\end{equation}
$\text{CRLB}(\hat{B})$는 데이터 증가에 $1/N^3$으로 감소하는 반면, $\text{CRLB}(\hat{A})$는 데이터 증가에 $1/N$ 비율로 감소한다. 이러한 차이로 인해 $x[n]$이 $B$를 변경하는 것에 $A$를 변경하는 것보다 더 민감하게 반응한다는 것을 알 수 있다.
3.5. Vector Parameter CRLB for Transformations
스칼라 파라미터의 변환에 대해 이전 섹션에서 알아본 것 처럼, 실제 추정 문제에서는 단순히 벡터 파라미터 $\boldsymbol{\theta}$를 추정해야 하는 것이 아닌 $\boldsymbol{\theta}$의 함수 형태를 추정해야 하는 일이 자주 발생한다.
추정하고자 하는 파라미터가 $\boldsymbol{\alpha} = \mathbf{g}(\boldsymbol{\theta})$와 같이 $\boldsymbol{\theta}$에 대한 함수이며 $r$ 차원의 함수의 함수일 경우 다음 수식을 만족해야 한다. (자세한 유도 과정은 Appendix 3B 참조)
\begin{equation} \label{eq:crlb9}
\begin{aligned}
\mathbf{C}_{\hat{\alpha}} - \frac{\partial \mathbf{g}(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} \mathbf{I}^{-1}(\boldsymbol{\theta}) \frac{\partial \mathbf{g}(\boldsymbol{\theta})^{\intercal}}{\partial \boldsymbol{\theta}} \geq 0
\end{aligned}
\end{equation}
- $(\cdot) \geq 0$ : 좌항에 있는 행렬이 positive semidefinite임을 의미한다.
자코비안 행렬 $\frac{\partial \mathbf{g}(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}}$는 $r\times p$ 크기의 행렬이며 다음과 같이 정의된다.
\begin{equation}
\begin{aligned}
\frac{\partial \mathbf{g}(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = \begin{bmatrix}
\frac{\partial g_{1}(\boldsymbol{\theta})}{\partial \theta_{1}} & \frac{\partial g_{1}(\boldsymbol{\theta})}{\partial \theta_{2}} & \cdots & \frac{\partial g_{1}(\boldsymbol{\theta})}{\partial \theta_{p}} \\
\frac{\partial g_{2}(\boldsymbol{\theta})}{\partial \theta_{1}} & \frac{\partial g_{2}(\boldsymbol{\theta})}{\partial \theta_{2}} & \cdots & \frac{\partial g_{2}(\boldsymbol{\theta})}{\partial \theta_{p}} \\
\vdots & \vdots & \ddots & \vdots \\
\frac{\partial g_{r}(\boldsymbol{\theta})}{\partial \theta_{1}} & \frac{\partial g_{r}(\boldsymbol{\theta})}{\partial \theta_{2}} & \cdots & \frac{\partial g_{r}(\boldsymbol{\theta})}{\partial \theta_{p}} \\
\end{bmatrix}
\end{aligned}
\end{equation}
3.5.1. Example 3.8 - CRLB for Signal-to-Noise Ratio
다음과 같은 관측 데이터가 주어졌다고 하자.
\begin{equation}
\begin{aligned}
x[n] = A + w[n] \quad \text{ where } w[n] \overset{\text{i.i.d.}}{\sim} \mathcal{N}(0, \sigma^{2})
\end{aligned}
\end{equation}
$\theta = [A, \ \sigma^2]^{\intercal}$는 미지의 파라미터이며 우리가 추정하고자 하는 값은 $\alpha = \frac{A^2}{\sigma^2}$라고 하자. 이러한 $\alpha$를 signal to noise ratio(SNR)이라고 한다. $\alpha = g(\boldsymbol{\theta}) = \theta_{1}^{2}/\theta_2 = A^2 / \sigma^2$과 같이 나타낼 수 있으므로 이를 통해 Fisher information 행렬을 표현하면 다음과 같다.
\begin{equation}
\begin{aligned}
\mathbf{I}(\boldsymbol{\theta}) = \begin{bmatrix}
\frac{N}{\sigma^{2}} & 0 \\ 0 & \frac{N}{2\sigma^{4}}
\end{bmatrix}
\end{aligned}
\end{equation}
자코비안 행렬 $\frac{\partial \mathbf{g}(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}}$은 다음과 같다.
\begin{equation}
\begin{aligned}
\frac{\partial \mathbf{g}(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} &= \begin{bmatrix}
\frac{\partial g(\boldsymbol{\theta})}{\partial \theta_{1}} & \frac{\partial g(\boldsymbol{\theta})}{\partial \theta_{2}} \\
\end{bmatrix} \\
&=\begin{bmatrix}
\frac{\partial g(\boldsymbol{\theta})}{\partial A} & \frac{\partial g(\boldsymbol{\theta})}{\partial \sigma^2} \\
\end{bmatrix} \\
& = \begin{bmatrix}
\frac{2A}{\sigma^{2}} & -\frac{A^2}{\sigma^4}
\end{bmatrix}
\end{aligned}
\end{equation}
따라서 (\ref{eq:crlb9}) 우측 항은 다음과 같이 구할 수 있다.
\begin{equation}
\begin{aligned}
\frac{\partial \mathbf{g}(\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} \mathbf{I}^{-1}(\boldsymbol{\theta}) \frac{\partial \mathbf{g}(\boldsymbol{\theta})^{\intercal}}{\partial \boldsymbol{\theta}} & = \begin{bmatrix}
\frac{2A}{\sigma^2} & -\frac{A^2}{\sigma^4}
\end{bmatrix}
\begin{bmatrix}
\frac{\sigma^2}{N} & 0 \\ 0 & \frac{2\sigma^2}{N}
\end{bmatrix}
\begin{bmatrix}
\frac{2A}{\sigma^2} \\ -\frac{A^2}{\sigma^4}
\end{bmatrix} \\
& = \frac{4A^2}{N\sigma^2} + \frac{2A^4}{N\sigma^4} \\
& = \frac{4\alpha + 2\alpha^2}{N}
\end{aligned}
\end{equation}
따라서 SNR 추정값 $\hat{\alpha}$에 대한 CRLB는 다음과 같다.
\begin{equation}
\begin{aligned}
\text{var}(\hat{\alpha}) \geq \frac{4\alpha + 2\alpha^2}{N}
\end{aligned}
\end{equation}
3.6. CRLB for the General Gaussian Case
이전 섹션에서 스칼라, 벡터 파라미터에 대한 CRLB를 알아보았다면 이번 섹션에서는 일반 가우시안 케이스에서 CRLB에 대해 알아본다. 일반 가우시안 케이스의 관측 데이터는 다음과 같이 나타낼 수 있다.
\begin{equation}
\begin{aligned}
\mathbf{x} \sim \mathcal{N}(\boldsymbol{\mu}(\boldsymbol{\theta}), \mathbf{C}(\boldsymbol{\theta}))
\end{aligned}
\end{equation}
위와 같이 평균과 표준편차는 파라미터 $\boldsymbol{\theta}$에 대해 종속적이다. 일반 가우시안 케이스에서 Fisher information 행렬은 다음과 같이 나타낼 수 있다. (자세한 유도 과정은 Appendix 3C 참조)
\begin{equation}
\begin{aligned}
\big[ \mathbf{I}(\boldsymbol{\theta}) \big]_{ij} & = \bigg[ \frac{\partial \boldsymbol{\mu}(\boldsymbol{\theta}) }{\partial \theta_{i}} \bigg]^{\intercal}
\mathbf{C}^{-1}(\boldsymbol{\theta})
\bigg[ \frac{\partial \boldsymbol{\mu}(\boldsymbol{\theta}) }{\partial \theta_{j}} \bigg] \\
& +\frac{1}{2} \text{tr} \bigg[ \mathbf{C}^{-1}(\boldsymbol{\theta}) \frac{\partial \mathbf{C}(\boldsymbol{\theta})}{\partial \theta_{i}} \mathbf{C}^{-1}(\boldsymbol{\theta}) \frac{\partial \mathbf{C}(\boldsymbol{\theta})}{\partial \theta_{j}} \bigg]
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
\text{where }\frac{\partial \boldsymbol{\mu}(\boldsymbol{\theta}) }{\partial \theta_{i}} = \begin{bmatrix}
\frac{\partial [\boldsymbol{\mu}(\boldsymbol{\theta})]_{1} }{\partial \theta_{i}} \\
\frac{\partial [\boldsymbol{\mu}(\boldsymbol{\theta})]_{2} }{\partial \theta_{i}} \\
\vdots \\
\frac{\partial [\boldsymbol{\mu}(\boldsymbol{\theta})]_{N} }{\partial \theta_{i}}
\end{bmatrix}
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
\text{ and }
\frac{\partial \mathbf{C}(\boldsymbol{\theta})}{\partial \theta_{i}} = \begin{bmatrix}
\frac{\partial [\mathbf{C}(\boldsymbol{\theta})]_{11} }{\partial \theta_{i}} & \frac{\partial [\mathbf{C}(\boldsymbol{\theta})]_{12} }{\partial \theta_{i}} & \cdots & \frac{\partial [\mathbf{C}(\boldsymbol{\theta})]_{1N} }{\partial \theta_{i}} \\
\frac{\partial [\mathbf{C}(\boldsymbol{\theta})]_{21} }{\partial \theta_{i}} & \frac{\partial [\mathbf{C}(\boldsymbol{\theta})]_{22} }{\partial \theta_{i}} & \cdots & \frac{\partial [\mathbf{C}(\boldsymbol{\theta})]_{2N} }{\partial \theta_{i}} \\
\vdots & \vdots & \ddots & \vdots \\
\frac{\partial [\mathbf{C}(\boldsymbol{\theta})]_{N1} }{\partial \theta_{i}} & \frac{\partial [\mathbf{C}(\boldsymbol{\theta})]_{N2} }{\partial \theta_{i}} & \cdots & \frac{\partial [\mathbf{C}(\boldsymbol{\theta})]_{NN} }{\partial \theta_{i}} \\
\end{bmatrix}
\end{aligned}
\end{equation}
스칼라 파라미터의 경우라면 관측 데이터는 다음과 같이 표현할 수 있다.
\begin{equation}
\begin{aligned}
\mathbf{x} \sim \mathcal{N}(\boldsymbol{\mu}(\theta), \mathbf{C}(\theta))
\end{aligned}
\end{equation}
Fisher information 행렬은 다음과 같이 간단하게 표현할 수 있다.
\begin{equation} \label{eq:crlb10}
\begin{aligned}
\mathbf{I}(\theta) & = \bigg[ \frac{\partial \boldsymbol{\mu}(\theta) }{\partial \theta} \bigg]^{\intercal}
\mathbf{C}^{-1}(\theta)
\bigg[ \frac{\partial \boldsymbol{\mu}(\theta) }{\partial \theta} \bigg] \\
& +\frac{1}{2} \text{tr} \bigg[ \mathbf{C}^{-1}(\theta) \frac{\partial \mathbf{C}(\theta)}{\partial \theta} \bigg]
\end{aligned}
\end{equation}
3.6.1. Example 3.11 - Random DC Level in WGN
다음과 같은 관측 데이터가 주어졌다고 가정하자.
\begin{equation}
\begin{aligned}
x[n] = A + w[n] \quad n=0,1,\cdots,N-1
\end{aligned}
\end{equation}
- $w[n] \sim \mathcal{N}(0, \sigma^2)$
- $A \sim \mathcal{N}(0, \sigma_{A}^2)$
관측 데이터 $\mathbf{x} = [x[0], x[1], \cdots, x[N-1]]^{\intercal}$이 주어졌을 때 이는 평균이 0이며 다음과 같은 $N\times N$ 크기의 공분산 행렬을 가진다.
\begin{equation}
\begin{aligned}
\big[ \mathbf{C}(\sigma_{A}^{2}) \big]_{ij} & = \mathbb{E}[x[i-1]x[j-1]] \\
& = \mathbb{E}[(A+w[i-1])(A+w[j-1])] \\
& = \sigma^{2}_{A} + \sigma^2 \delta_{ij}
\end{aligned}
\end{equation}
\begin{equation}
\begin{aligned}
\therefore \mathbf{C}(\sigma_{A}^{2}) & = \sigma^{2}_{A} \mathbf{11}^{\intercal} + \sigma^2 \mathbf{I}
\end{aligned}
\end{equation}
- $\mathbf{1} = [1, 1, \cdots, 1]^{\intercal}$
Woodbury identity 공식에 의해 공분산의 역행렬은 다음과 같다. 자세한 내용은 해당 링크를 참조하면 된다.
\begin{equation}
\begin{aligned}
\mathbf{C}^{-1}(\sigma_{A}^{2}) & = \frac{1}{\sigma^2} \bigg( \mathbf{I} - \frac{\sigma^2_{A}}{\sigma^2 + N\sigma_{A}^{2}} \mathbf{11}^{\intercal} \bigg)
\end{aligned}
\end{equation}
공분산의 편미분은 다음과 같다.
\begin{equation}
\begin{aligned}
\frac{ \partial \mathbf{C}(\sigma_{A}^{2}) }{\partial \sigma_A^2} = \mathbf{11}^{\intercal}
\end{aligned}
\end{equation}
따라서 둘을 합치면 다음 공식이 성립한다.
\begin{equation}
\begin{aligned}
\mathbf{C}^{-1}(\sigma_{A}^{2}) \frac{ \partial \mathbf{C}(\sigma_{A}^{2}) }{\partial \sigma_A^2} = \frac{1}{\sigma^2 + N \sigma_A^2} \mathbf{11}^{\intercal}
\end{aligned}
\end{equation}
최종적으로 (\ref{eq:crlb10}) 식은 다음과 같이 쓸 수 있다.
\begin{equation}
\begin{aligned}
\mathbf{I}(\sigma_{A}^{2}) & = \frac{1}{2} \text{tr} \bigg[
\bigg( \frac{1}{\sigma^2 + N \sigma^{2}_{A}} \bigg)^2 \mathbf{11}^{\intercal} \mathbf{11}^{\intercal} \bigg] \\
& = \frac{N}{2} \bigg( \frac{1}{\sigma^2 + N \sigma^{2}_{A}} \bigg)^2 \text{tr}(\mathbf{11}^{\intercal}) \\
& = \frac{1}{2} \bigg( \frac{N}{\sigma^2 + N \sigma^{2}_{A}} \bigg)^2
\end{aligned}
\end{equation}
따라서 $\sigma_{A}^{2}$에 대한 CRLB는 다음과 같이 나타낼 수 있다.
\begin{equation}
\begin{aligned}
\text{var}(\sigma_{A}^{2}) \geq 2 \bigg( \sigma_{A}^{2} + \frac{\sigma^2}{N} \bigg)^{2}
\end{aligned}
\end{equation}
위 식에서 보다시피 $N \rightarrow \infty$가 되어도 CRLB는 $2\sigma_{A}^{4}$ 이하로 내려가지 않는다. 이는 각각의 개별 관측 데이터가 전부 $A$에 대한 정보를 포함하고 있기 때문에 $A$의 분산 값이 항상 반영되기 때문이다.
4. Linear Models
일반적으로 MVUE를 찾는 과정은 쉬운 작업이 아니다. 하지만 다행이도 많은 신호 처리 문제들은 선형(linear) 데이터 모델의 형태를 가지는데 이러한 특성이 쉽게 MVUE를 결정할 수 있도록 도와준다. 선형 모델임이 확인되면 MVUE 값이 즉시 명백해지는 것 뿐만 아니라 통계적인 성능도 자연스럽게 따라온다. 그러므로 최적의 추정값을 찾는 열쇠는 문제를 선형 모델로 구조화하여 그 고유한 특성을 활용하는 것이 중요하다.
4.1. Definition and Properties
다음과 같은 line fitting 문제가 주어졌다고 하자.
\begin{equation}
\begin{aligned}
x[n] = A + Bn + w[n] \quad n=0,1,\cdots,N-1
\end{aligned}
\end{equation}
- $w[n] \sim \mathcal{N}(0, \sigma^{2})$ : WGN
이 때, y절편의 값 $A$와 직선의 기울기 $B$의 값을 찾아야 한다. 위 식을 행렬 형태로 작성하면 다음과 같다.
\begin{equation} \label{eq:lm1}
\begin{aligned}
\mathbf{x} = \mathbf{H}\boldsymbol{\theta} + \mathbf{w}
\end{aligned}
\end{equation}
각 원소를 펼쳐보면 다음과 같다.
\begin{equation}
\begin{aligned}
& \mathbf{x} = [x[0],x[1],\cdots,x[N-1]]^{\intercal} \\
& \mathbf{w} = [w[0],w[1], \cdots, w[N-1]]^{\intercal} \\
& \boldsymbol{\theta} = [A, \ B ]^{\intercal} \\
& \mathbf{H} = \begin{bmatrix}
1&0 \\ 1 & 1 \\ \vdots & \vdots \\ 1 & N-1
\end{bmatrix}
\end{aligned}
\end{equation}
행렬 $\mathbf{H}$는 $N\times2$의 크기를 가지며 일반적으로 관측 함수(observation matrix)라고 부른다. 데이터 $\mathbf{x}$는 파라미터 $\boldsymbol{\theta}$가 관측 행렬 $\mathbf{H}$를 통과한 다음에 관측되기 때문이다. 노이즈 벡터도 $\mathbf{w} \sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I})$와 같이 벡터 형태로 나타낼 수 있다. 이러한 (\ref{eq:lm1})의 형태를 선형 모델(liinear model)이라고 부른다. 선형모델의 노이즈는 일반적으로 가우시안 형태를 가진다.
앞서 3장의 Theorm 3.2에서 본 것처럼 $\hat{\boldsymbol{\theta}} = \mathbf{g}(\mathbf{x})$가 MVUE가 되려면 다음 식을 만족해야 한다.
\begin{equation} \label{eq:lm2}
\begin{aligned}
\frac{\partial \ln p(\mathbf{x};\boldsymbol{\theta}) }{\partial \boldsymbol{\theta}} = \mathbf{I}(\boldsymbol{\theta})(\mathbf{g}(\mathbf{x}) - \boldsymbol{\theta})
\end{aligned}
\end{equation}
그리고 위 식을 만족할 때 $\hat{\boldsymbol{\theta}}$의 분산은 $\mathbf{I}^{-1}(\boldsymbol{\theta})$가 된다. 이러한 Theorem 3.2의 조건들을 선형 모델에 대해 확장하면 다음과 같다.
\begin{equation}
\begin{aligned}
\frac{\partial \ln p(\mathbf{x};\boldsymbol{\theta}) }{\partial \boldsymbol{\theta}} & = \frac{\partial }{\partial \boldsymbol{\theta}} \bigg[
-\ln(2\pi\sigma^2)^{\frac{N}{2}} - \frac{1}{2\sigma^2}(\mathbf{x} - \mathbf{H}\boldsymbol{\theta})^{\intercal}(\mathbf{x} - \mathbf{H}\boldsymbol{\theta}) \bigg] \\
& = -\frac{1}{2\sigma^2} \frac{\partial }{\partial \boldsymbol{\theta}} [\mathbf{x}^{\intercal}\mathbf{x} - 2\mathbf{x}^{\intercal}\mathbf{H}\boldsymbol{\theta} + \boldsymbol{\theta}^{\intercal}\mathbf{H}^{\intercal}\mathbf{H}\boldsymbol{\theta}]
\end{aligned}
\end{equation}
위 식에 다음과 같은 등식을 적용한다.
\begin{equation}
\begin{aligned}
& \frac{\partial \mathbf{b}^{\intercal}\boldsymbol{\theta}}{\partial \boldsymbol{\theta}} = \mathbf{b} \\
& \frac{\partial \boldsymbol{\theta}^{\intercal}\mathbf{A}\boldsymbol{\theta}}{\partial \boldsymbol{\theta}} = 2\mathbf{A}\boldsymbol{\theta}
\end{aligned}
\end{equation}
- $\mathbf{A}$ : 임의의 대칭(symmetric) 행렬
\begin{equation}
\begin{aligned}
\frac{\partial \ln p(\mathbf{x};\boldsymbol{\theta}) }{\partial \boldsymbol{\theta}} = \frac{1}{\sigma^2} [\mathbf{H}^{\intercal}\mathbf{x} - \mathbf{H}^{\intercal}\mathbf{H}\boldsymbol{\theta} ]
\end{aligned}
\end{equation}
만약 $\mathbf{H}^{\intercal}\mathbf{H}$의 역행렬이 존재한다면 식은 다음과 같다.
\begin{equation}
\begin{aligned}
\frac{\partial \ln p(\mathbf{x};\boldsymbol{\theta}) }{\partial \boldsymbol{\theta}} = \frac{\mathbf{H}^{\intercal}\mathbf{H}}{\sigma^2} [(\mathbf{H}^{\intercal}\mathbf{H})^{-1}\mathbf{H}^{\intercal}\mathbf{x} - \boldsymbol{\theta} ]
\end{aligned}
\end{equation}
위 식은 (\ref{eq:lm2})의 형태와 정확히 일치한다.
\begin{equation}
\boxed{ \begin{aligned}
& \hat{\boldsymbol{\theta}} = (\mathbf{H}^{\intercal}\mathbf{H})^{-1}\mathbf{H}^{\intercal}\mathbf{x} \\
& \mathbf{I}(\boldsymbol{\theta}) = \frac{\mathbf{H}^{\intercal}\mathbf{H}}{\sigma^2}
\end{aligned} }
\end{equation}
따라서 MVUE를 만족하는 $\hat{\boldsymbol{\theta}}$의 분산은 $\mathbf{I}^{-1}(\boldsymbol{\theta})$가 된다.
\begin{equation}
\boxed{ \begin{aligned}
& \mathbf{C}_{\hat{\boldsymbol{\theta}}} = \mathbf{I}^{-1}(\boldsymbol{\theta}) = \sigma^2 (\mathbf{H}^{\intercal}\mathbf{H})^{-1}
\end{aligned} }
\end{equation}
\begin{equation}
\boxed{ \begin{aligned}
\therefore \hat{\boldsymbol{\theta}} \sim \mathcal{N}(\boldsymbol{\theta}, \sigma^2 (\mathbf{H}^{\intercal}\mathbf{H})^{-1})
\end{aligned} }
\end{equation}
4.1.1. Example 4.2 - Fourier Analysis
많은 신호는 주기적인 특성을 지닌다. 이러한 강한 주기적 특성은 푸리에 해석(Fourier analysis)를 통해 수학적으로 모델링할 수 있다. 푸리에 계수(coefficient)록 해당 주파수 성분이 강하게 포함되어 있다는 의미를 지닌다. 해당 예제에서는 푸리에 해석을 하는 과정이 선형 모델 파라미터를 추정하는 것과 동일한 과정임을 보인다. 다음과 같은 정현파 신호에 가우시안 노이즈가 추가된 예제를 보자.
\begin{equation}
\begin{aligned}
x[n] = \sum_{k=1}^{M} a_k \cos\bigg( \frac{2\pi k n}{N}\bigg) + \sum_{k=1}^{M} b_k \sin\bigg( \frac{2\pi k n}{N}\bigg) + w[n] \quad n = 0,1,\cdots, N-1
\end{aligned}
\end{equation}
- $w[n] \sim \mathcal{N}(0, \sigma^{2})$ : WGN
주파수 성분은 $f_1 = 1/N$부터 $f_k =k/N$ 성분이 복합적으로 섞여 있는 것으로 가정한다. 여기서 우리는 $a_k, b_k$ 계수를 추정해야 한다. 벡터 파라미터를 다음과 같이 정의한다.
\begin{equation}
\begin{aligned}
\boldsymbol{\theta} = [a_1, a_2, \cdots, a_M, b_1, b_2, \cdots, b_M]^{\intercal}
\end{aligned}
\end{equation}
그리고 $\mathbf{H}$ 행렬을 다음과 같이 정의한다.
\begin{equation}
\begin{aligned}
\mathbf{H} = \begin{bmatrix}
1 & \cdots & 1 & 0 & \cdots & 0 \\
\cos(\frac{2\pi}{N}) & \cdots & \cos(\frac{2\pi M}{N}) & \sin(\frac{2\pi}{N}) & \cdots & \sin(\frac{2\pi M}{N}) \\
\vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\
\cos(\frac{2\pi(N-1)}{N}) & \cdots & \cos(\frac{2\pi M(N-1)}{N}) & \sin(\frac{2\pi(N-1)}{N}) & \cdots & \sin(\frac{2\pi M(N-1)}{N})
\end{bmatrix}
\end{aligned}
\end{equation}
위 식에서 행렬의 크기는 $N \times 2M$이며 $p=2M$이 된다. $\mathbf{H}$가 $N > p$를 만족하려면 $M < N/2$가 되어야 한다. $\mathbf{H}$의 각 열벡터는 DFT 특성에 의해 서로 직교(orthogonal)하므로 계산을 편하게 하기 위해 이를 열벡터로 표시한다.
\begin{equation}
\begin{aligned}
\mathbf{H} = \begin{bmatrix}
\mathbf{h}_{1} & \mathbf{h}_{2} & \cdots & \mathbf{h}_{2M}
\end{bmatrix}
\end{aligned}
\end{equation}
직교성(orthogonality)에 의해 다음 식을 만족한다.
\begin{equation}
\begin{aligned}
\mathbf{h}_{i}^{\intercal} \mathbf{h}_{j} =0 \quad \text{ for } i \neq j
\end{aligned}
\end{equation}
$\mathbf{H}^{\intercal}\mathbf{H}$를 전개해보면 다음과 같다.
\begin{equation}
\begin{aligned}
\mathbf{H}^{\intercal}\mathbf{H} & = \begin{bmatrix}
\mathbf{h}_{1}^{\intercal} \\ \vdots \\ \mathbf{h}_{2M}^{\intercal}
\end{bmatrix}
\begin{bmatrix}
\mathbf{h}_{1} & \cdots & \mathbf{h}_{2M}
\end{bmatrix} \\
& = \begin{bmatrix}
\mathbf{h}_{1}^{\intercal}\mathbf{h}_{1} & \mathbf{h}_{1}^{\intercal}\mathbf{h}_{2} & \cdots & \mathbf{h}_{1}^{\intercal}\mathbf{h}_{2M} \\
\mathbf{h}_{2}^{\intercal}\mathbf{h}_{1} & \mathbf{h}_{2}^{\intercal}\mathbf{h}_{2} & \cdots & \mathbf{h}_{2}^{\intercal}\mathbf{h}_{2M} \\
\vdots & \vdots & \ddots & \vdots \\
\mathbf{h}_{2M}^{\intercal}\mathbf{h}_{1} & \mathbf{h}_{2M}^{\intercal}\mathbf{h}_{2} & \cdots & \mathbf{h}_{2M}^{\intercal}\mathbf{h}_{2M} \\
\end{bmatrix}
\end{aligned}
\end{equation}
위 식은 직교성 특징으로 인해 대각 행렬(diagonal matrix)가 되고 매우 쉽게 역행렬을 구할 수 있다. DFT의 직교 성질은 다음과 같다.
\begin{equation}
\begin{aligned}
& \sum_{n=0}^{N-1} \cos \bigg( \frac{2\pi i n}{N} \bigg) \cos \bigg( \frac{2\pi j n}{N} \bigg) = \frac{N}{2} \delta_{ij} \\
& \sum_{n=0}^{N-1} \sin \bigg( \frac{2\pi i n}{N} \bigg) \sin \bigg( \frac{2\pi j n}{N} \bigg) = \frac{N}{2} \delta_{ij} \\
& \sum_{n=0}^{N-1} \cos \bigg( \frac{2\pi i n}{N} \bigg) \sin \bigg( \frac{2\pi j n}{N} \bigg) = 0 \quad \text{ for all } i,j
\end{aligned}
\end{equation}
직교성 특징을 적용하면 다음과 같이 간단한 식으로 변한다.
\begin{equation}
\begin{aligned}
\mathbf{H}^{\intercal}\mathbf{H} & = \begin{bmatrix}
\frac{N}{2} & 0 & \cdots & 0 \\
0 & \frac{N}{2} & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & \frac{N}{2} \\
\end{bmatrix} = \frac{N}{2} \mathbf{I}
\end{aligned}
\end{equation}
$\hat{\boldsymbol{\theta}}$가 MVUE이면 다음 식을 만족해야 한다.
\begin{equation}
\begin{aligned}
\hat{\boldsymbol{\theta}} & = (\mathbf{H}^{\intercal}\mathbf{H})^{-1}\mathbf{H}^{\intercal}\mathbf{x} \\
& = \frac{2}{N} \mathbf{H}^{\intercal}\mathbf{x} \\
& = \frac{2}{N} \begin{bmatrix}\mathbf{h}_{1}^{\intercal} \\ \vdots \\ \mathbf{h}_{2M}^{\intercal} \end{bmatrix} \mathbf{x} \\
& = \begin{bmatrix} \frac{2}{N}\mathbf{h}_{1}^{\intercal}\mathbf{x} \\ \vdots \\ \frac{2}{N}\mathbf{h}_{2M}^{\intercal}\mathbf{x} \end{bmatrix}
\end{aligned}
\end{equation}
따라서 $a_k, b_k$는 다음과 같이 추정할 수 있다.
\begin{equation}
\boxed{ \begin{aligned}
&\hat{a_{k}} = \frac{2}{N} \sum_{n=0}^{N-1} x[n] \cos \bigg( \frac{2\pi k n}{N} \bigg) \\
& \hat{b_{k}} = \frac{2}{N} \sum_{n=0}^{N-1} x[n] \sin \bigg( \frac{2\pi k n}{N} \bigg)
\end{aligned} }
\end{equation}
이는 DFT 계수와 동일하다. 추정값의 기대값과 분산을 구해보면 다음과 같다.
\begin{equation}
\boxed{ \begin{aligned}
\mathbb{E}(\hat{a_{k}}) & = a_k \\
\mathbb{E}(\hat{b_{k}}) & = b_k \\
\mathbf{C}_{\hat{\boldsymbol{\theta}}} & = \sigma^2 (\mathbf{H}^{\intercal}\mathbf{H})^{-1} \\
& = \sigma^2 \big( \frac{N}{2} \mathbf{I} \big)^{-1} \\
& = \frac{2\sigma^2}{N} \mathbf{I}
\end{aligned} }
\end{equation}
추정값의 분산 $\mathbf{C}_{\hat{\boldsymbol{\theta}}}$이 대각 행렬이므로 각 파라미터들은 서로 독립적이다. 이러한 독립적 특성으로 인해 $\mathbf{H}$의 열벡터가 서로 직교하여 $\mathbf{H}^{\intercal}\mathbf{H}$가 매우 단순하게 계산되었다. 만약 주파수 성분을 임의의 성분으로 변경하면 직교성이 상실되어 MVUE를 구하는 과정이 매우 복잡해진다.
4.2. Extension to the Linear Model
일반적인 경우 선형 모델의 노이즈는 WGN이 아니다. 선형 모델을 일반적인 경우로 확장하면 노이즈는 다음과 같이 모델링할 수 있다.
\begin{equation}
\begin{aligned}
\mathbf{w} \sim \mathcal{N}(0, \mathbf{C})
\end{aligned}
\end{equation}
$\mathbf{C}$ 행렬은 더 이상 스칼라 값이 곱해진 대각 identity 행렬이 아니다. 노이즈는 0 이상의 값을 가지기 때문에 $\mathbf{C}$는 positive definite 행렬이라고 가정할 수 있고 따라서 $\mathbf{C}^{-1}$ 또한 positive definite이다. 이를 다음과 같이 분해할 수 있다.
\begin{equation}
\begin{aligned}
\mathbf{C}^{-1} = \mathbf{D}^{\intercal}\mathbf{D}
\end{aligned}
\end{equation}
$\mathbf{D}$ 행렬은 $N \times N$ 크기의 역행렬이 존재하는 행렬이다. $\mathbf{D}$ 행렬은 기존 노이즈를 whitening 변환해주는데 사용된다.
\begin{equation}
\begin{aligned}
\mathbb{E}[(\mathbf{D}\mathbf{w})(\mathbf{D}\mathbf{w})^{\intercal}] & = \mathbf{DCD}^{\intercal} \\
& = \mathbf{D}\mathbf{D}^{-1}\mathbf{D}^{-\intercal}\mathbf{D}^{\intercal} \\
& = \mathbf{I}
\end{aligned}
\end{equation}
일반적인 선형 모델은 다음과 같다.
\begin{equation}
\begin{aligned}
\mathbf{x} = \mathbf{H}\boldsymbol{\theta} + \mathbf{w}
\end{aligned}
\end{equation}
여기에 whitening을 위해 다음과 같은 변환을 거친다.
\begin{equation}
\begin{aligned}
\mathbf{x}' & = \mathbf{Dx} \\
& = \mathbf{DH} \boldsymbol{\theta} + \mathbf{Dw} \\
& = \mathbf{H}' \boldsymbol{\theta} + \mathbf{w}'
\end{aligned}
\end{equation}
위 식은 $\mathbf{w}' \sim \mathcal{N}(0, \mathbf{I})$와 같이 whitening 되었다. 다음 순서는 앞서 설명한 선형 모델 공식과 동일하다.
\begin{equation}
\begin{aligned}
\hat{\boldsymbol{\theta}} & = (\mathbf{H}'^{\intercal}\mathbf{H}')^{-1}\mathbf{H}'^{\intercal}\mathbf{x}' \\
& = (\mathbf{H}^{\intercal}\mathbf{D}^{\intercal}\mathbf{D}\mathbf{H})^{-1}\mathbf{H}^{\intercal}\mathbf{D}^{\intercal}\mathbf{D}\mathbf{x} \\
\end{aligned}
\end{equation}
\begin{equation}
\boxed{ \begin{aligned}
\hat{\boldsymbol{\theta}} & = (\mathbf{H}^{\intercal}\mathbf{C}^{-1}\mathbf{H})^{-1}\mathbf{H}^{\intercal}\mathbf{C}^{-1}\mathbf{x}
\end{aligned} }
\end{equation}
분산도 동일한 과정을 거친다.
\begin{equation}
\begin{aligned}
\mathbf{C}_{\hat{\boldsymbol{\theta}}} & = (\mathbf{H}'^{\intercal}\mathbf{H}')^{-1}
\end{aligned}
\end{equation}
\begin{equation}
\boxed{ \begin{aligned}
\mathbf{C}_{\hat{\boldsymbol{\theta}}} & = (\mathbf{H}'^{\intercal}\mathbf{C}^{-1}\mathbf{H}')^{-1}
\end{aligned} }
\end{equation}
Wrap up
References