9주차-5월 02일 (1)
빅데이터분석특강
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
$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
- suppose: $p=0.2$ 이라고 하자.
그렇다면 $(x_1,x_2,x_3,x_4)=(0,1,0,1)$와 같은 샘플이 얻어질 확률이 아래와 같다.
0.8 * 0.2 * 0.8 * 0.2
- 질문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])
- 확률이라는 말을 쓸 수 없지만 확률의 느낌은 있음 -> 가능도라는 말을 쓰자.
- 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}$
- 이 예제의 경우 $p$의 최대우도추정량을 구하면
$$\hat{p}^{mle} = \text{argmax}_p L(p) = \text{argmax}_p \big\{ p^2(1-p)^2 \big\}= \frac{1}{2}$$
(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
- 문제를 풀면서 생각해보니 손실함수는 -로그가능도함수로 선택하면 될 것 같다?
- 손실함수를 선택하는 기준이 -로그가능도함수만 존재하는 것은 아니나 대부분 그러하긴함
(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
(해석) - 로지스틱에서 가능도함수와 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
- 회귀분석이든 로지스틱이든 손실함수는 minus_log_likelihood 로 선택한다.
- 그런데 (오차항이 정규분포인) 회귀분석 일때는 minus_log_likelihood 가 MSEloss가 되고
- 로지스틱일때는 minus_log_likelihood 가 BCEloss가 된다
- minus_log_likelihood가 손실함수를 선택하는 유일한 기준은 아니다. <--- 참고만하세요, 이 수업에서는 안중요합니다.
- 오차항이 대칭이고 서로독립이며 등분산 가정을 만족하는 어떠한 분포에서의 회귀모형이 있다고 하자. 이 회귀모형에서 $\hat{\beta}$은 여전히 MSEloss를 최소화하는 $\beta$를 구함으로써 얻을 수 있다.
- 이 경우 MSEloss를 쓰는 이론적근거? $\hat{\beta}$이 BLUE가 되기 때문임 (가우스-마코프정리)