최대 우도 추정 소개 (in R)

데이터 분석/확률 통계

최대 우도 추정 소개 (in R)

yessen 2023. 1. 3. 19:00

728x90

# 본 포스팅은 Analytics Vidhya의 An Introductory Guide to Maximum Likelihood Estimation (with a case study in R)을 번역한 것입니다.

최대 우도 추정 (Maximum Likelihood Expectation)

예를 들어 회사의 주가를 예측하는 모델을 구축했다고 가정 해 보자.

밤새 주가가 급격히 상승한 것을 관찰했다면, 그 뒤에는 여러 가지 이유가 있을 수 있다.

가장 가능성이 높은 이유의 가능성을 찾는 것이 Maximum Likelihood Estimation의 전부이다.

이 게시물에서는 최대 우도 추정(이하 MLE)이 작동하는 방식과 분포를 사용하여 모델의 계수를 결정하는 데 사용할 수있는 방법을 살펴 본다.

참고 : 이 내용은 수학의 기초와 확률을 알고 있다고 가정함.

Why should I use Maximum Likelihood Estimation (MLE)?

이벤트 티켓 판매를 예측한다고 가정 해 보자. 데이터는 다음과 같은 히스토그램과 밀도를 가짐.

변수를 어떻게 모델링하면 좋을까? 여기에서 변수는 정규 분포가 아니며 비대칭이므로 선형 회귀 가정을 위반한다.

널리 사용되는 방법은 log, sqrt, reciprocal 등으로 변수를 변환하여,

변환된 변수가 정규 분포를 따르도록 하여 선형 회귀로 모델링하는 것이다.

각 변환을 시도해서 결과를 살펴보자.

With Log transformation:

With Square Root Transformation:

With Reciprocal:

이들 중 어느 것도 정규 분포에 가깝지 않다. (--> ~~정규화 가정 실패~~)

모델의 기본 가정을 위반하지 않도록 하려면 이러한 데이터를 어떻게 모델링해야 할까?

정규 분포가 아닌 다른 분포로 모델링하는 것은 어떨까? 다른 분포를 사용한다면 계수는 어떻게 추정할까?

MLE(Maximum Likelihood Estimation)는 이러한 경우 잇점을 제공한다.

Understanding MLE with an example

통계와 확률을 공부하다보면 다음과 같이 부딪치는 문제가 있다.

x가 평균 50과 표준 편차 (sd) 10을 갖는 정규 분포를 따르는 경우 x > 100의 확률은 얼마일까?

우리는 이미 분포(이 경우 정규 분포) 및 해당 파라미터(평균 및 sd)를 알고 있지만, 실제 문제에서는 이러한 정보를 알 수 없으며 데이터에서 추정해야 한다.

MLE는 주어진 데이터를 가장 잘 설명하는 분포의 모수를 결정하는 데 도움이되는 기술임.

예를 들어, 학생들의 체중 (kg)을 나타내는 데이터가 있다고 가정하자.

x = as.data.frame(rnorm(50,50,10))
ggplot(x, aes(x = x)) + geom_dotplot()

이것은 정규 분포를 따르는 것으로 보이지만, 이 분포에 대한 평균 및 표준 편차 (sd)는 어떻게 얻을 수 있을까?

한 가지 방법은 주어진 데이터의 평균과 sd를 직접 계산하는 것이며, 이는 각각 49.8Kg와 11.37이다.

이 값은 주어진 데이터를 잘 나타내지만 모집단을 가장 잘 설명하지는 못할 수도 있다.

이때, 보다 강력한 모수 추정값을 얻기 위해 MLE을 사용할 수 있다.

즉, MLE는 관찰된 데이터를 획득할 확률(가능성)이 최대화되도록 샘플 데이터로부터 모집단 파라미터 (예를 들어, 평균에 대한 평균 및 분산, 포아송에 대한 속도 (lambda) 등)를 추정하는 방법을 위해 정의될 수 있다.

MLE를 직관적으로 이해하려면, 다음 중 위 그림의 데이터를 관찰할 확률을 최대화하는 것을 추측해 보아라.

평균 = 100, SD = 10

평균 = 50, SD = 10

확실히 모집단 평균이 100인 경우 위의 데이터 형태를 관찰할 가능성은 거의 없다.

Getting to know the technical details

MLE가 할 수있는 것에 대한 직관을 얻었으므로 실제로 가능성과 그 최대화 방법에 대해 자세히 알아보도록 하자.

먼저 분포 모수에 대한 빠른 검토부터 시작하겠다.

Distribution Parameters

먼저 분포 모수를 이해해 보자. 이 용어에 대한 Wikipedia의 정의는 “확률 분포의 계열을 가리키는 수량”이다.

즉, 모집단의 수치적 특성 또는 통계 모델로 간주할 수 있으며, 다음 다이어그램으로 이해할 수 있다.

벨 커브의 너비와 높이는 평균과 분산이라는 두 가지 매개 변수에 의해 결정되는데, 이것을 정규 분포의 분포 모수라고 한다. 마찬가지로, 포아송 분포는 람다 (lambda)라는 하나의 매개 변수에 의해 제어되는데, 람다 (lambda)는 이벤트가 시간 또는 공간 간격으로 발생하는 횟수를 의미한다.

대부분의 분포에는 하나 또는 두 개의 파라미터가 있지만 일부 분포는 4 파라미터 베타 분포와 같이 최대 4개의 파라미터가 있을 수도 있다.

Likelihood

위 그림에서 분포 파라미터 세트가 주어지면 일부 데이터 값이 다른 데이터보다 가능성이 높다는 것을 알 수 있다.

맨 처음 그림에서, 우리는 주어진 데이터가 평균이 100이 아닌 50 일 때 발생할 가능성이 더 높다는 것을 보았지만, 이것은 사실 우리가 이미 데이터를 관찰했기에 알 수 있었다.

따라서 우리는 역의 문제에 직면하게 된다 -> 관측된 데이터와 관심있는 모델이 주어졌을 때, 우리는 이러한 데이터를 생산할 확률이 높은 모든 확률 밀도 중에서 하나의 확률 밀도 함수/확률 질량 함수 (f (x | θ))를 찾아야 한다.

이 역 문제를 해결하기 위해 데이터 벡터 x의 역할과 f (x | θ)의 (분포) 모수 벡터 θ의 역할을 반대로하여 우도 함수를 정의한다. 즉,

L (θ; x) = f (x | θ)

MLE에서, 우도 함수 L(θ; x)이 있다고 가정 할 수 있는데, 여기서 θ는 분포 모수 벡터이고 x는 관측값 세트이다.

우리는 주어진 관측치 (x의 값)로 가능성을 최대화하는 θ 값을 찾는 데 관심이 있다.

Log Likelihood

관측치 (xi)가 확률 분포 f0 (여기서 f0 = 정규 분포)에서 도출된 독립적이고 동일하게 분포된 랜덤 변수라고 가정하면 수학적 문제는 더 단순해 진다. 이렇게 하면 우도 함수가 다음과 같이 줄어든다.

이 함수의 최대값 / 최소값을 찾기 위해 이 함수의 미분 값을 w.r.t θ로 취하여 0에 대한 방정식으로 풀 수 있다(기울기 0은 최대값 또는 최소값을 나타냄).

여기에 내적(product)에 대한 항이 있으므로 내적에 체인 규칙(chain rule)을 적용해야 한다.

이를 좀 더 쉽게 하기 위해 우도 함수에 로그를 취하고, 똑같이 이를 최대화 하한다. (로그는 곱을 합으로 변환할 수 있고, 로그는 증가하는 함수이므로 결과 값 θ에 영향을 미치지 않기 때문)

그러므로 다음을 구할 수 있다.

Maximizing the Likelihood

로그 우도 함수 LL (θ; x)의 최대 값을 찾기 위해 다음을 수행할 수 있다.

1. LL (θ; x) 함수 w.r.t θ의 1 차 미분 값을 취하여 0에 대한 방정식을 푼다

2. LL (θ; x) 함수 w.r.t θ의 2 차 미분 값을 취하여 음의 값인지 확인한다

미적분이 우도를 최대화하는 데 직접적인 도움이 되지 않는 상황도 많이 있지만, 최대 값은 쉽게 확인할 수 있다.

로그 가능성을 최대화하는 매개 변수 값을 찾을 때 어떤 종류의 '우선 순위'나 특별한 위치에 1 차 도함수를 0으로 설정하는 것은 없고, 단지 몇 가지 매개 변수를 추정해야 할 때 편리한 도구이다.

일반적으로, 함수의 argmax를 식별하기 위한 거의 모든 유효한 접근법이 로그 우도 함수의 최대값을 찾는 데 적합 할 수 있는데, 이것은 제한되지 않은 비선형 최적화 문제이다. 우리는 다음과 같은 방식으로 동작하는 최적화 알고리즘을 찾게 된다.

임의의 시작점에서 local minia로 수렴
가능한 빨리 실행되는 것

우도를 최대화하기 위해서는 최적화 기술을 사용하는 것이 일반적이며, 다양한 방법이 있다 (뉴턴 방법, 피셔 스코어링, 다양한 켤레 기울기 기반 접근 방법, 가장 가파른 하강, Nelder-Mead 유형 (단순) 접근 방법, BFGS 및 기타 다양한 기술).

위의 예에서와 같이 모델이 가우스 인 것으로 가정하면 MLE 추정치는 일반적인 최소 제곱법(OLS)과 같다.

Determining model coefficients using MLE

이제 MLE을 사용하여 예측 모델의 계수를 결정하는 방법을 살펴 보자.

Yi ~ P (µi)를 사용하여 독립적인 포아송 랜덤 변수 n개의 관측치 y1, y2 . . , yn의 표본이 있다고 가정해 보자.

또, 평균 µi (따라서 분산!)가 설명 변수 xi의 벡터에 의존한다고 가정하면, 다음과 같이 간단한 선형 모델을 만들 수 있다.

여기서 θ는 모델 계수의 벡터이다. 이 모형은 우항의 선형 예측 변수는 어떤 실제 값이라도 가정할 수 있는 반면, 예상 횟수를 나타내는 좌항의 포아송 평균은 음수가 아니어야한다는 단점이 있다. 이 문제에 대한 간단한 해결책은 선형 모델을 사용하여 평균의 로그를 모델링하는 것이다. 따라서 로그로 연결된 일반화 선형 모델을 다음과 같이 만들 수 있다.

우리의 목표는 MLE을 사용하여 θ를 찾는 것.

이제 포아송 분포는 다음과 같다.

θ를 찾기 위해 이전 섹션에서 배운 로그 우도 개념을 적용 할 수 있다.

위 방정식에 로그를 취하고 log (y!)와 관련된 상수를 무시하면 로그 우도 함수는 다음과 같다.

여기서 μi는 공변량 xi와 θ 계수에 의존한다. µi = exp (xi'θ)를 대입하여 방정식을 풀면 우도를 최대화하는 θ를 구할 수 있고, θ 벡터가 있으면 xi와 θ 벡터를 곱하여 평균(기대)값을 예측할 수 있다.

MLE using R

이 섹션에서는 실제 데이터 셋을 사용하여 앞에서 배운 개념을 사용하여 문제를 해결해 보겠다. (데이터 다운로드 링크)

Datetime Count of tickets sold

25-08-2012 00:00 8

25-08-2012 01:00 2

25-08-2012 02:00 6

25-08-2012 03:00 2

25-08-2012 04:00 2

25-08-2012 05:00 2

이 데이터는 2012 년 8 월 25 일부터 2014 년 9 월 25 일까지 (약 18K 레코드) 매시간 판매된 티켓 수를 보유하고 있다. 우리의 목표는 각 시간에 판매된 티켓 수를 예측하는 것이다.

이 문제는 회귀, 시계열 등과 같은 기술을 사용하여 해결할 수 있는데,

여기서는 R을 사용하여 위에서 배운 통계 모델링 기술을 사용한다.

먼저, 통계 모델링에서 우리는 목표(target) 변수가 어떻게 분포되어 있는지에 관심이 있다.

카운트 분포를 살펴 보자.

hist(Y$Count, breaks = 50,probability = T ,main = "Histogram of Count Variable")
lines(density(Y$Count), col="red", lwd=2)

이것은 포아송 분포나 지수 분포로도 피팅 할 수 있어 보인다.

현재 변수는 티켓 수이므로 Poisson이 더 적합한 모델이다

(지수 분포는 일반적으로 이벤트 간의 시간 간격을 모델링하는 데 사용됨).

2년 동안 판매된 티켓 수를 플랏팅하면 다음과 같다.

시간이 지남에 따라 티켓 판매가 크게 증가한 것이 보이며, 일을 단순하게 하기 위해 연령만을 요인으로 하여 결과를 모델링하면 다음과 같다. 여기서 연령은 2012년 8월 25일 이후의 주(week) 수이다.

여기서 µ (판매된 티켓 수)는 포아송 분포의 평균을 따르는 것으로 가정하고 θ0 및 θ1은 추정해야하는 계수이다.

Eq1, Eq2를 결합하면 다음과 같이 로그 우도 함수를 얻을 수 있다.

R stats4 패키지에서 mle() 함수를 사용하여 계수 θ0 및 θ1을 추정 할 수 있으며, 다음과 같은 기본 파라미터가 필요하다.

최소화해야 할 음수 가능도 함수 : 이것은 방금 도출 한 것과 같지만 앞에 음수 부호가 있다 (로그 가능성을 최대화하는 것은 음수 로그 가능성을 최소화하는 것과 같다)
계수 벡터의 시작점 : 계수의 초기 추측 값. (함수가 로컬 최소값에 도달 할 수 있으므로 결과는 이러한 값에 따라 달라질 수 있으며, 다른 시작점으로도 실행하여 결과를 확인하는 것이 좋다)
선택적으로, 최적화 되어야 하는 우도 함수를 사용하는 방법. BFGS가 기본 방법임.

이 예에서는 음수 로그 우도 함수를 다음과 같이 코딩 할 수 있다.

nll <- function(theta0,theta1) {    
	x <- Y$age[-idx]    
	y <- Y$Count[-idx]    
    mu = exp(theta0 + x*theta1) - sum(y*(log(mu)) - mu)
    }

모델의 성능을 객관적으로 평가할 수 있도록 데이터를 train / test로 나누었다.

idx는 테스트 세트에있는 행의 인덱스.

set.seed(200)
idx <- createDataPartition(Y$Count, p=0.25,list=FALSE)

다음으로 mle 함수를 호출하여 파라미터를 가져옴.

est <- stats4::mle(minuslog=nll, start=list(theta0=2,theta1=0))
summary(est)

Maximum likelihood estimation
Call:stats4::mle(minuslogl = nll, start = list(theta0 = 2, theta1 = 0))
Coefficients:         
		Estimate  Std. Error
theta0  2.68280754  2.548367e-03
theta1  0.03264451  2.998218e-05

-2 log L: -16594396

이로써 계수의 추정치를 알 수 있음.

테스트 세트에 대한 결과를 얻기 위해 RMSE를 평가 지표로 사용.

pred.ts <- (exp(coef(est)['theta0'] + Y$age[idx]*coef(est)['theta1'] ))
rmse(pred.ts, Y$Count[idx])

86.95227

이제 이 모델을 카운트 로그로 모델링 된 표준 선형 모델 (오류가 정규 분포인)과 비교해 보겠다.

lm.fit <-  lm(log(Count)~age, data=Y[-idx,])
Coefficients:             
			Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.9112992 0.0110972   172.2 <2e-16 ***age         
			0.0414107 0.0001768   234.3 <2e-16 ***
            
pred.lm <- predict(lm.fit, Y[idx,])
rmse(exp(pred.lm), Y$Count[idx])

93.77393

표준 선형 모형에 대한 RMSE는 포아송 분포가있는 모형보다 높다.

두 모델의 잔차 플랏을 비교해, 어떻게 모델의 성능이 다른 영역에서 차이가 나는지 확인해 보자.

정규 선형 회귀와 비교했을 때, 포아송 회귀를 사용한 에러가 0에 훨신 더 가까운 것을 볼 수 있다.

Python에서는 BFGS, L-BFGS등의 방법과 파라미터의 최소화나 초기값 추정을 위해 scipy.optimize.minimize()을 사용하면 이와 같은 결과를 얻을 수 있다.

R에서는 stat 패키지에서 glm을 사용하여 집단 분포를 더 간단히 모델링할 수 있으며, Poisson, Gamma, Binomial, Quasi, Inverse Gaussian, Quasi Binomial, Quasi Poisson 분포를 모두 지원한다.

위의 예제와 같이, 아래와 같은 명령어를 사용하면 계수를 직접적으로 얻을 수 있다.

glm(Count ~ age, family = "poisson", data = Y) 

Coefficients:             
	Estimate Std. Error z value Pr(>|z|)     
(Intercept) 2.669   2.218e-03    1203 <2e-16 *** 
age 		0.03278 2.612e-05    1255 <2e-16 ***

저작자표시 (새창열림)