9주차-5월 02일 (1)

imports

import numpy as np
import tensorflow as tf
import tensorflow.experimental.numpy as tnp

tnp.experimental_enable_numpy_behavior()

import matplotlib.pyplot as plt

우도함수와 최대우도추정량

예제

$X_i \overset{iid}{\sim} Ber(p)$에서 얻은 샘플이 아래와 같다고 하자.

x=[0,1,0,1]
x

[0, 1, 0, 1]

$p$는 얼마라고 볼 수 있는가? --> 0.5

왜?? $p$가 0.5라고 주장할 수 있는 이론적 근거, 혹은 논리체계가 무엇인가?

- suppose: $p=0.1$ 이라고 하자.

그렇다면 $(x_1,x_2,x_3,x_4)=(0,1,0,1)$와 같은 샘플이 얻어질 확률이 아래와 같다.

0.9 * 0.1 * 0.9 * 0.1

0.008100000000000001

- suppose: $p=0.2$ 이라고 하자.

그렇다면 $(x_1,x_2,x_3,x_4)=(0,1,0,1)$와 같은 샘플이 얻어질 확률이 아래와 같다.

0.8 * 0.2 * 0.8 * 0.2

0.025600000000000008

- 질문1: $p=0.1$인것 같냐? 아니면 $p=0.2$인것 같냐? -> $p=0.2$

• 왜?? $p=0.2$일 확률이 더 크다!

(여기서 잠깐 중요한것) 확률이라는 말을 함부로 쓸 수 없다.

- 0.0256은 "$p=0.2$일 경우 샘플 (0,1,0,1)이 얻어질 확률"이지 "$p=0.2$일 확률"은 아니다.

"$p=0.2$인 확률" 이라는 개념이 성립하려면 아래코드에서 sum([(1-p)*p*(1-p)*p for p in _plist])이 1보다는 작아야 한다. (그런데 1보다 크다)

_plist = np.linspace(0.499,0.501,1000)
sum([(1-p)*p*(1-p)*p for p in _plist])

62.49983299986714

- 확률이라는 말을 쓸 수 없지만 확률의 느낌은 있음 -> 가능도라는 말을 쓰자.

• 0.0256 $=$ $p$가 0.2일 경우 샘플 (0,1,0,1)이 얻어질 확률 $=$ $p$가 0.2일 가능도

- 다시 질문1로 돌아가자!

• 질문1: $p=0.1$인 것 같냐? 아니면 $p=0.2$인 것 같냐? -> 답 $p=0.2$ -> 왜? $p=0.2$인 가능도가 더 크니까!
• 질문2: $p=0.2$인 것 같냐? 아니면 $p=0.3$인 것 같냐? -> 답 $p=0.3$ -> 왜? $p=0.3$인 가능도가 더 크니까!
• 질문3: ...

- 궁극의 질문: $p$가 뭐일 것 같아?

• $p$가 입력으로 들어가면 가능도가 계산되는 함수를 만들자.
• 그 함수를 최대화하는 $p$를 찾자.
• 그 $p$가 궁극의 질문에 대한 대답이 된다.

- 잠깐 용어정리

• 가능도함수 $=$ 우도함수 $=$ likelihood function $:=$ $L(p)$
• $p$의 maximum likelihood estimator $=$ p의 MLE $:=$ $\hat{p}^{mle}$ $=$ $\text{argmax}_p L(p)$ $=$ $\hat{p}$

(예제의 풀이)

- 이 예제의 경우 가능도함수를 정의하자.

• $L(p)$: $p$의 가능도함수 = $p$가 모수일때 샘플 (0,1,0,1)이 얻어질 확률 = $p$가 모수일때 $x_1$이 0일 확률 $\times \dots \times$ $p$가 모수일때 $x_4$가 1일 확률
• $L(p)=\prod_{i=1}^{4} f(x_i;p)= \prod_{i=1}^{4}p^{x_i}(1-p)^{1-x_i}$

Note: 참고로 이 과정을 일반화 하면 $X_1,\dots,X_n \overset{iid}{\sim} Ber(p)$ 일때 $p$의 likelihood function은 $\prod_{i=1}^{n}p^{x_i}(1-p)^{1-x_i}$ 라고 볼 수 있다.

Note: 더 일반화: $x_1,\dots,x_n$이 pdf가 $f(x)$인 분포에서 뽑힌 서로 독립인 샘플일때 likelihood function은 $\prod_{i=1}^{n}f(x_i)$라고 볼 수 있다.

- 이 예제의 경우 $p$의 최대우도추정량을 구하면

$$\hat{p}^{mle} = \text{argmax}_p L(p) = \text{argmax}_p \big\{ p^2(1-p)^2 \big\}= \frac{1}{2}$$

중간고사 1번

(1) $N(\mu,\sigma)$에서 얻은 샘플이 아래와 같다고 할때 $\mu,\sigma$의 MLE를 구하여라.

<tf.Tensor: shape=(10000,), dtype=float64, numpy=
array([ 4.12539849,  5.46696729,  5.27243374, ...,  2.89712332,
        5.01072291, -1.13050477])>

(2) $Ber(p)$에서 얻은 샘플이 아래와 같다고 할 때 $p$의 MLE를 구하여라.

<tf.Tensor: shape=(10000,), dtype=int64, numpy=array([1, 1, 1, ..., 0, 0, 1])>

(3) $y_i = \beta_0 + \beta_1 x_i + \epsilon_i$, $\epsilon_i \overset{iid}{\sim} N(0,1)$ 일때 $(\beta_0,\beta_1)$의 MLE를 구하여라. (회귀모형)

(풀이) 가능도함수

$$L(\beta_0,\beta_1)=\prod_{i=1}^{n}f(y_i), \quad f(y_i)=\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}(y_i-\mu_i)^2}, \quad \mu_i=\beta_0+\beta_1 x_i$$

를 최대화하는 $\beta_0,\beta_1$을 구하면된다. 그런데 이것은 아래를 최소화하는 $\beta_0,\beta_1$을 구하는 것과 같다.

$$-\log L(\beta_0,\beta_1) = \sum_{i=1}^{n}(y_i-\beta_0-\beta_1x_i)^2$$

위의 식은 SSE와 같다. 결국 오차항이 정규분포를 따르는 회귀모형의 MLE는 MSE를 최소화하는 $\beta_0,\beta_1$을 구하면 된다.

중간고사 1-(3)의 다른 풀이

step1: 생성

x= tf.constant(np.arange(1,10001)/10000)
y= tnp.random.randn(10000) + (0.5 + 2*x)

step2: minimize MSEloss (원래는 maximize log-likelihood)

• maximize likelihood였던 문제를 minimize MSEloss로 바꾸어도 되는근거? 주어진 함수(=가능도함수)를 최대화하는 $\beta_0,\beta_1$은 MSE를 최소화하는 $\beta_0,\beta_1$과 동일하므로

beta0= tf.Variable(1.0)
beta1= tf.Variable(1.0)
for i in range(2000):
    with tf.GradientTape() as tape:
        #minus_log_likelihood = tf.reduce_sum((y-beta0-beta1*x)**2)
        loss =  tf.reduce_sum((y-beta0-beta1*x)**2)
    slope1, slope2 = tape.gradient(loss,[beta0,beta1])
    beta0.assign_sub(slope1* 0.1/10000) # N=10000
    beta1.assign_sub(slope2* 0.1/10000)

beta0,beta1

(<tf.Variable 'Variable:0' shape=() dtype=float32, numpy=0.47993016>,
 <tf.Variable 'Variable:0' shape=() dtype=float32, numpy=2.020918>)

- 문제를 풀면서 생각해보니 손실함수는 -로그가능도함수로 선택하면 될 것 같다?

• 손실함수를 선택하는 기준이 -로그가능도함수만 존재하는 것은 아니나 대부분 그러하긴함

(4) 출제하지 못한 중간고사 문제

아래의 모형을 생각하자.

• $Y_i \overset{iid}{\sim} Ber(\pi_i)$
• $\pi_i = \frac{\exp(w_0+w_1x_i)}{1+\exp(w_0+w_1x_i)}=\frac{\exp(-1+5x_i)}{1+\exp(-1+5x_i)}$

아래는 위의 모형에서 얻은 샘플이다.

x = tnp.linspace(-1,1,2000)
pi = tnp.exp(-1+5*x) / (1+tnp.exp(-1+5*x))
y = np.random.binomial(1,pi)
y = tf.constant(y)

함수 $L(w_0,w_1)$을 최대화하는 $(w_0,w_1)$를 tf.GradeintTape()를 활용하여 추정하라. (경사하강법 혹은 경사상승법을 사용하고 $(w_0,w_1)$의 초기값은 모두 0.1로 설정할 것)

$$L(w_0,w_1)=\prod_{i=1}^{n}f(y_i), \quad f(x_i)={\pi_i}^{y_i}(1-\pi_i)^{1-y_i},\quad \pi_i=\text{sigmoid}(w_0+w_1x_i)$$

(풀이1)

w0hat = tf.Variable(1.0)
w1hat = tf.Variable(1.0)

for i in range(1000):
    with tf.GradientTape() as tape:
        pihat = tnp.exp(w0hat+w1hat *x) / (1+tnp.exp(w0hat+w1hat *x))
        pdf = pihat**y * (1-pihat)**(1-y)
        logL = tf.reduce_mean(tnp.log(pdf))
    slope1,slope2 = tape.gradient(logL,[w0hat,w1hat])
    w0hat.assign_add(slope1*0.1)
    w1hat.assign_add(slope2*0.1)

w0hat,w1hat

(<tf.Variable 'Variable:0' shape=() dtype=float32, numpy=-0.88931483>,
 <tf.Variable 'Variable:0' shape=() dtype=float32, numpy=4.231925>)

(해석) - 로지스틱에서 가능도함수와 BCEloss의 관계

$L(w_0,w_1)$를 최대화하는 $w_0,w_1$은 아래를 최소화하는 $w_0,w_1$와 같다.

$$-\log L(w_0,w_1) = - \sum_{i=1}^{n}\big(y_i \log(\pi_i) + (1-y_i)\log(1-\pi_i)\big)$$

이것은 최적의 $w_0,w_1$을 $\hat{w}_0,\hat{w}_1$이라고 하면 $\hat{\pi}_i=\frac{\exp(\hat{w}_0+\hat{w}_1x_i)}{1+\exp(\hat{w}_0+\hat{w}_1x_i)}=\hat{y}_i$이 되고 따라서 위의 식은 $n\times$BCEloss의 형태임을 쉽게 알 수 있다.

결국 로지스틱 모형에서 $(w_0,w_1)$의 MLE를 구하기 위해서는 BCEloss를 최소화하는 $(w_0,w_1)$을 구하면 된다!

(풀이2)

w0hat = tf.Variable(1.0)
w1hat = tf.Variable(1.0)

for i in range(1000):
    with tf.GradientTape() as tape:
        yhat = tnp.exp(w0hat+w1hat *x) / (1+tnp.exp(w0hat+w1hat *x))
        loss = tf.losses.binary_crossentropy(y,yhat)
    slope1,slope2 = tape.gradient(loss,[w0hat,w1hat])
    w0hat.assign_sub(slope1*0.1)
    w1hat.assign_sub(slope2*0.1)

w0hat,w1hat

(<tf.Variable 'Variable:0' shape=() dtype=float32, numpy=-0.88931495>,
 <tf.Variable 'Variable:0' shape=() dtype=float32, numpy=4.2319255>)

손실함수의 설계 (선택)

- 회귀분석이든 로지스틱이든 손실함수는 minus_log_likelihood 로 선택한다.

• 그런데 (오차항이 정규분포인) 회귀분석 일때는 minus_log_likelihood 가 MSEloss가 되고
• 로지스틱일때는 minus_log_likelihood 가 BCEloss가 된다

- minus_log_likelihood가 손실함수를 선택하는 유일한 기준은 아니다. <--- 참고만하세요, 이 수업에서는 안중요합니다.

• 오차항이 대칭이고 서로독립이며 등분산 가정을 만족하는 어떠한 분포에서의 회귀모형이 있다고 하자. 이 회귀모형에서 $\hat{\beta}$은 여전히 MSEloss를 최소화하는 $\beta$를 구함으로써 얻을 수 있다.
• 이 경우 MSEloss를 쓰는 이론적근거? $\hat{\beta}$이 BLUE가 되기 때문임 (가우스-마코프정리)