ABOUT ME

-

Today
-
Yesterday
-
Total
-
  • 7장 최대 우도 추정법
    ML DL DS/사회과학자를위한데이터과학 2022. 3. 16. 12:08

    7장 최대 우도 추정법

    로날드 피셔님이 최대우도추정법을 완성함.

    MLE - 수리통계학 (참고)


    • Likelihood Function 이란

      세타로 부터 추출된, 독립인 표본들의 결합분포로 보는 것이다.

      x가 주어질 때, 세타에 대한 함수로 본다.

    • Regularity Condition - 사전 규칙

      R0: 세타가 다르면, pdf도 다르다

      R1: pdf들은 공통의 support를 가진다. pdf가 support란 0보다 큰 값을 갖게 되는 영역

      ex. 균등분포면 해당 구간, 정규분포면 -\infty \sim \infty

      R2: 세타의 참값 세타 제로는 오메가에 있는 interior point이다.

      세타 값의 open set을 잡아도 오메가에 들어간다.

      interior point: 경계가 아닌 점 [a, b]에서는 a b를 제외한 a< x < b 상의 점들

    • Likelihood의 성질
      • 실제 파라미터: 표본의 크기가 충분히 클 때, θ0\theta_{0}의 가능도가 다른 θ\theta의 가능도보다 클 확률을 θ0\theta_{0}에서 evaluate할 때, 1로 가까이 간다.
        • 증명

          첫 째줄의 우변에서 로그 안을 어떤 변수로 보면 기댓값 → Weak Law of Large Numbers

          E[log(z)] < log E[z] → 옌센의 부등식, cf) https://engineer-mole.tistory.com/90

          기댓값 계산에 의해 세타 0의 pdf를 곱해주면 분모만 남아 pdf의 기댓값은 1 → log 1 = 0

    • MLE란

      보통 아래로 구함, 아니면 미분을 통해 작아지는 방향으로 품.

      (항상 맞으려면, 두 번 미분해서 음수이면, 최대라는 것을 알 수 있음. 파라미터가 p 개일 경우: Negative Definite → concave하다는 것의 정의)

      Negative Definite: https://blog.naver.com/yjh0853/221676307215

    • MLE 예시
      • (i) (iI) 지수, 라플라스분포: explicit solution이 있는 경우
      • (iii) 로지스틱의 경우: explicit solution은 없지만, MLE가 있는지 확인하고 이를 numerical 하게 구할 수 있는지를 확인. 미분한 식을 정리해서 k(세타)를 만들고, 미분하면 항상 양수이기 때문에 k(세타)는 항상 증가하는 것을 알 수 있다. 세타가 -무한대로 갈 때 k(세타)는 0으로, 세타가 무한대로 갈때는 n으로 수렴한다. MLE는 K(세타) = 2/n 가 되게 하는 세타가 MLE이다. 이를 grid search를 통해 구할 수 있다.
      • (iv) (v) 균등분포, 베르누이분포: 제약조건이 있어야 구할 수 있는 경우
        • 균등분포: https://www.youtube.com/watch?v=clEVsBgu-f4

          균등분포의 모수가 0~세타가지라는 제약 조건

          균등분포의 가능도함수 그래프를 보면, 세타의 추정치로 표본 중 최대로 나온 표본 값으로 설정할 때 크다

        • 베르누이분포:

          베르누이분포의 모수가 0보다크고 1/3보다 작다는 제약조건

    • MLE 예시 2
    • MLE 의 성질
      • Functional Invariance

        MLE는 일대일 함수에 대하여 변하지 않는 성질

        함수 g 가 일대일 함수이고, MLE of θ=θ^\theta = \hat{\theta} 이면 → MLE of g(θ)=g(θ^)g(\theta) = g(\hat{\theta}) 이다

        ex. θ=Xˉ\theta = \bar{X}eθ=eXˉe^{\theta} = e^{\bar{X}}, log(θ)=log(Xˉ)\log({\theta}) = \log({\bar{X}}) ....

      • Consistency

        로그우도함수의 세타에 대한 편미분함수 = 0에 대한 해에 대한 추정치 세타 햇은 n이 커질 수록, 실제 모수에 확률수렴(converge in probability) 한다.

    최대우도추정법을 이용한 통계적 추론


    • 로그우도를 많이 사용하는 이유: 대부분의 확률 분포가 log 변형 시 concave 해짐
    • 최대우도추정법은 베이지안 방법과 마찬가지로 확률의 전도(inversion of probability)관점에서 통계적 추론을 접근
      • 확률의 전도: 관측된 자료를 통해 알지 못하는 값을 학습
      • 통계적 추론: 자료를 통해 자료의 생성과정을 학습, 추론 대상은 자료의 생성과정을 설명하는 모수
      • MLE의 속성

        일관성: 관측자료가 무한으로 증가하면 확률적으로 추정값은 모수와 가까워진다.

        근사적 정규성(aympototic normality): 관측자료가 무한으로 증가하면 모수값으로 1/n1/\sqrt{n} 의 속도로 수렴

      • 두 속성은 MLE가 베이지안 방법과 대등한 통계적 추론 방법이면서, 베이지안 처럼 사전분포를 필요로 하지 않는다는 점에 주목을 받음
      • 통계적 추론에는 기대값과 기대값의 불확실성에 대한 정보가 필요

        불확실성 정보: 서베이 - 표본오차, 빈도주의 통계 - 신뢰구간, 베이지안 방법 - 확률구간(Bayesian credible interval)

        MLE - 피셔정보한도 (Fisher Information Bound, I(θ)I(\theta))로 불확실성을 측정함

      • MLE의 불확실성: Fisher Information, score function
        • Score Function: 로그우도함수의 1차도함수:

          s(θ)=(logL(θ;x))θ=θlogf(x;θ)=f(x;θ)f(x;θ)s(\theta) = {{\partial{(\log{L(\theta;x)})}}\over{\partial \theta}} = {{\partial}\over{\partial \theta}} \log{f(x;\theta)} = {{f'(x;\theta)}\over{f(x;\theta)}}

        • (Expected) Fisher Information: score function의 분산

          I(θ;x)=V[θlogf(x;θ)]I(\theta;x) = V[{{\partial}\over{\partial \theta}} \log{f(x;\theta)}] 

          =((logf(x;θ))θ)2=E[((logf(X;θ))θ)2]= \int_{-\infty}^{\infty}({{\partial{(\log{f(x;\theta)})}}\over{\partial \theta}})^{2} = E[({{\partial{(\log{f(X;\theta)})}}\over{\partial \theta}})^{2}]

          =(2logf(x;θ)θ2)=E[(2logf(X;θ)θ2)]= - \int_{-\infty}^{\infty} ({{\partial^{2}\log{f(x;\theta)} }\over{\partial\theta^{2}}}) = - E[({{\partial^{2}\log{f(X;\theta)} }\over{\partial\theta^{2}}})]

          • 의미
            • 표본마다 가능도는 아래 그림과 같이 다르다.
            • 각 표본마다의 가능도의 기울기가 score 함수
            • 참모수에 대한 최대우도추정치 세타스타 과 표본들에 대한 스코어함수의 기댓값 → 0
            • 피셔정보량의 의미: 기울기들의 분산.

              기울기의 분산이 크다 → 그래프가 뾰족함이 크다는 것을 의미함 → 정보가 많음, 표본수가 많아지면 분산(피셔정보량)이 커짐

            • Observed Fisher Information은 참모수 대신에 최대우도추정치를 기준으로 표본들의 가능도 기울기 분산을 계산

            ref. https://www.youtube.com/watch?v=wjsNqBmjBuw

            • 모수는 모집단의 고정값이지만, 스코어함수의 기댓값이 0이라고 가정하는 경우는, 모집단의 모수가 확률변수인 경우에 논리적 일관성을 갖게될 수 있음

    일반 선형모형


    • 가우스-마르코프 정리를 이용한 최소제곱 추정치는 오차가 정규분포를 따르는 것을 가정
    • 하지만 사회과학 자료는 항상 실수 형태를 취하지 않음
      • 종속변수의 사건형 자료형, 양의 실수값만 취하는 자료(시간 등), 범주형 자료 등
    • 이러한 종속변수에 대해 선형회귀모형을 적용할 경우 문제 발생 가능성이 있음
      1. 가우스-마르코프 가정 위반 → 점추정치 편향 발생, 구간 추정치 불확실성 정보가 점추정치의 불확실성을 정확하게 표현하지 못함
      1. 추정치의 해석이 어려워짐 → 독립변수의 단위 증가에 대한 종속 변수의 영향
      1. 해석의 어려움으로 인한 예측의 어려움
    • 일반 선형 모형 등장(McCullagh and Nelder(1989))
      • 핵심 아이디어: 제한된 종속변수를 실수형태로 변환(변수변환)해 선형 모형에 연결
      • 단순하지만 효과가 좋음
      • 이분형 / 다분형 자료(dischotomous / polychotomous), 분절/단절(truncated/censored), 시간지속 자료(duration data)등을 종속변수로 갖는 회귀모형에 적용 가능
    • 일반 선형 모형을 구성하는 3가지 요소
      • 확률분포

        종속변수의 특성을 고려해 자료 생성과정을 설명하는 확률분포

        일반 선형 모형의 확률분포는 지수분포족(정규분포, 이항분포, 포아송분포 등)에서 선택됨

        지수분포족의 확률분포함수는 아래로 표현 가능. ϕ\phi는 산포도를 결정하는 산포모수(dispersion parameter), f는 밀도함수:

        f(y;θ,ϕ)=exp[θyb(θ)a(ϕ)+c(y,ϕ)]f(y;\theta,\phi) =\exp[{{\theta y - b(\theta)}\over{a(\phi)}}+c(y,\phi)]

        기댓값: E(Yi)=b(θi)E(Y_{i}) = b'(\theta_{i})

        분산: V(Yi)=a(ϕ)b(θi)V(Y_{i}) = a(\phi)b''(\theta_{i})

      • 선형함수

        설명변수와 파라미터의 선형결합(XβX\beta)은 종속변수의 경향성을 설명하는 은닉변수

      • 연결함수 (일반선형모형을 구성하는 핵심요소)

        확률분포의 기댓값(μ\mu)을 선형함수에 연결하는 함수 g

        g(μ)=Xβg(\mu) = X\beta

        평균함수: 연결함수의 역함수, 선형함수를 기댓값에 연결하는 역할

        μ=g1(Xβ)\mu = g^{-1}(X\beta)

    • 요소들간의 관계

      추정은 오른쪽에서 왼쪽으로 가는 과정

      추정 목표 값: β^\hat\beta , 예측과 진단은 베타햇을 통해 Y를 다시 생성하는 것. Y_true와 Y_pred를 이용해 모형의 타당성 확인가능.

    최대우도추정법을 이용한 일반 선형 모형의 추정


    일반 선형모형 r의 glm은 반복가중최소제곱법(iteratively reweighted least squares)를 사용해 최대우도추정값을 찾음

    θj+1=arg minθΩwi(θj)Σi=1nyifi(θ)p\theta^{j+1} = \argmin_{\theta \in \Omega} w_{i}(\theta^{j})\Sigma_{i=1}^{n}|y_{i}-f_{i}(\theta)|^{p}

    • 푸아송분포의 그리드 탐색을 이용한 MLE
      • 푸아송 분포: 단위 시간 안에 어떤 사건이 몇 번 발생할 것인지를 표현하는 이산확률분포

        (람다: 정해진 시간 안에 어떤 사건이 일어날 횟수에 대한 기댓값, k: 그 사건이 k회 일어날 확률)

      • 푸아송분포의 가능도 함수

        L(λ;y)=i=1nexp(λ)λyiyi!L(\lambda ;y) = \prod_{i=1}^{n} {{\exp(-\lambda) \lambda^{y_{i}}}\over{y_{i}!}}

      • 푸아송분포의 로그 가능도 함수

        l(λ;y)=Σi=1n(yilog(λ)λlog(yi!))l(\lambda ;y) = \Sigma_{i=1}^{n} (y_{i}\log(\lambda) - \lambda - \log({y_{i}!}))

      • 평균을 3으로 갖는 푸아송분포에서 10000개 표본 추출, 1~6의 100개의 값들에 대하여 로그가능도를 확인해 람다 추정
        • R에서 그리드 탐색을 통한 MLE 예시
      • 그리드 탐색은 파라미터의 수가 적고, 탐색 범위가 제한적이다. 추정해야하는 파라미터가 늘어나고 경사하강이 일관되지 않은 모습을 보이는 경우 뉴튼 랩슨과 같은 방법을 사용하는 것이 일반적인 MLE 추정법이다.

    • 뉴튼 랩슨 방법을 이용한 MLE ( newton-raphson method)
      • 참고- 뉴튼 랩슨 방법

        현재 x값에서 접선을 그리고 접선이 x축과 만나는 지점으로 x를 이동시켜 가면서 점진적으로 해를 찾는 방법. f(x)=0의 해를 수치적으로 찾을 때 유용 → f'(x)=0을 이용해 가능도 함수 f(x)의 최대값을 찾자

        기하학적으로: xt+1x_{t+1} 를 기울기의 역수에 높이를 곱한 값으로 xtx_{t}을 이동 시키자

        zero-finding 문제 → 원래는 zero-finding 문제에 해결하는 것인데, 최적화의 필요조건인 f’(x)=0인 지점을 찾을 때도 활용

        딥러닝의 손실함수를 최적화하는 방법인 Gradient Descent 방법을 대체가능하고, learning rate 부분이 없어 더 정확하고 빠르게 계산할 수 있다. 하지만, 파라미터가 너무 많아 2차 헤시안까지 구해야하는 계산 비용 문제, 오버피팅문제 등이 있을 수 있다. 이를 위해 일차 미분만으로 계산이 가능한 가우스 뉴튼법을 사용할 수 있다. (연립인 경우 가우스 뉴튼, 쿼사이 뉴턴은 gradient의 평균을 이용해 hessian을 근사하는 알고리즘)

        xt+1=xtf(xt)f(xt)x_{t+1} = x_{t} - {{f(x_{t})}\over{f'(x_{t})}}

        뉴튼 랩슨을 MLE에서 최대 문제 적용할 경우, (f'(x)=0을 이용해 가능도 함수 f(x)의 최대값을 찾기)에 이용할 때에는 아래와 같이 변형됨

        xt+1=xtf(xt)f(xt)x_{t+1} = x_{t} - {{f'(x_{t})}\over{f''(x_{t})}}

      • 선형 회귀분석 모형의 최대우도추정

        추정해야할 파라미터가 많아질 수록 그리드 탐색은 어렵다

        α=1,β=1,σϵ2=1\alpha =1, \beta=-1, \sigma_{\epsilon}^{2}=1을 모수로 갖는 단순 선형 회귀분석 모형을 통해 Y를 생성하고, X는 균등분포에서 생성해, α,β,σϵ2\alpha, \beta, \sigma_{\epsilon}^{2} (y절편, 기울기, 잔차의 분산)에 대해 MLE를 수행해보자.

        lm.like <- function(theta, y, X){ 
          n <- nrow(X)
          k <- ncol(X)
          beta <- theta[1:k] 
          sigma <- exp(theta[k + 1])
          loglike <- sum(log(dnorm(y, X %*% beta, sigma)))	## dnorm: 정규분포의 pdf 
          return(loglike)
        }
        # ⍺=1, β=-1, σ^2=1 인 데이터 생성
        set.seed(1999)
        X <- cbind(1, runif(100)) 
        theta.true <- c(1,-1, 1) # 정답 데이터 생성
        y <- X%*%theta.true[1:2] + rnorm(100)
        
        plot(x =X[,2], y= y)
        
        mle <- optim(par=c(5,5,5), fn=lm.like, method="BFGS", control = list(fnscale = -1), 
                   hessian=TRUE, y=y, X=X)
        
        ## 표준편차 
        K <- ncol(X)
        psi.mu <- mle$par[K + 1]
        # solve(행렬) -> 역행렬
        psi.sd <- sqrt(solve(-mle$hessian)[K + 1, K + 1])
        sigma.mle <- exp(psi.mu + 0.5*psi.sd) 
        sigma.mle # 0.9143931
        
        ols <- lm(y ~ X-1)
        summary(ols)$sigma # 0.89159
        • optim 함수

          control = list(fnscale = -1): 최댓값을 찾으라는 것

          hessian=True: logf(y;θ)\log f''(y;\theta)를 보고될 객체에 포함하라는 것

          par=c(5,5,5): 파라미터 초기값

        • 추정 결과

          α=\alpha =1.15, β=\beta= -1.1로 추정

          이 때, 표준편차의 경우 로그정규분포변환 때문에 계산이 간단하지 않다.

          ZN(μ,σ2)Z \sim \mathcal{N}(\mu,\sigma^{2})의 지수화 할 때(Ui=exp(Zi)U_{i} = \exp(Z_{i}))의 평균은 아래와 같다.

          E(Ui)=exp(μ+σ/2)E(U_{i}) = \exp(\mu + \sigma/2 )

          cf. https://ko.wikipedia.org/wiki/로그_정규_분포

          MLE를 통해 구한 잔차의 표준편차 0.9143은 ols를 통해 구한 잔차의 표준편차 0.8915와 유사함을 알수 있다.

        import numpy as np
        from scipy.stats import norm
        from scipy.optimize import minimize
        import pandas as pd
        import statsmodels.formula.api as smf
        
        # R로 생성한 데이터 
        data = pd.read_csv('./data.csv')
        X = data[['X1','X2']]
        X = np.array(X)
        theta_true = [1, -1, 1]
        
        
        def likelihood(theta,y,X):
            n = len(X)
            k = len(X[0])
            beta = theta[0:k]
            sigma = np.exp(theta[k])
            log_like = np.sum(np.log(norm(X @ beta, sigma).pdf(y)))
            return -log_like
        
        mle = minimize(fun=likelihood, 
                 x0=np.array([5,5,5]), 
                 args=(y, X),
                 # hess=False,
                 method='BFGS', 
                 options={'disp': True})
        # x: array([ 1.15352384, -1.10766175, -0.12485006])
        
        print(np.exp(mle.x[2] + 0.5 * mle.hess_inv[2][2]))
        # 0.8849
        
        model = smf.ols('y ~ X1+X2 -1 ', data =data)
        result=model.fit()
        print(result.summary())
        print()
        print(np.sqrt(result.scale)) # 0.89159 
        print(result.resid.std(ddof=2)) # 0.89159

      • 프로빗 모형 MLE - R example

        프로빗 모형의 확률분포는 아래와 같다. 종속변수가 {0,1}의 이분변수일 때 사용한다.

        조건부 분포를 YXBin(y;μ(x;β),N)Y|X \sim Bin(y;\mu(x;\beta),N)라고 설정하며

        μ(x)\mu(x)를 log odds함수의 역함수인 로지스틱함수로 평균함수σ()\sigma()로 설정하면 로지스틱회귀이고, 표준정규분포의 cdf인 Φ()\Phi()로 평균함수로 설정하면 probit model 이라 부른다.

        • 한 샘플에 대한 확률: (x’: x transpose를 의미)
        • 로그 우도
        • gradient, Lβ{{\partial\mathcal{L}}\over{\partial\beta}}

          =Σi=1N(ϕ(xiTβ)(yΦ(xiTβ))Φ(xiTβ)(1Φ(xiTβ)))= \Sigma_{i=1}^{N} ({{\phi(x^{T}_{i} \beta) (y-\Phi(x^{T}_{i}\beta))}\over{\Phi(x^{T}_{i}\beta)(1- \Phi(x^{T}_{i}\beta))}})

          , where ϕ\phi: cdf of standard normal distribution

        • 데이터

          low :indicator of birth weight less than 2.5 kg.

          age : mother's age in years.

          lwt : mother's weight in pounds at last menstrual period.

          race : mother's race (1 = white, 2 = black, 3 = other).

          smoke: smoking status during pregnancy.

          ptl: number of previous premature labours.

          ht: history of hypertension.

          ui: presence of uterine irritability.

          ftv: number of physician visits during the first trimester.

          bwt: birth weight in grams.

        suppressMessages(library(MCMCpack))
        data(birthwt)
        
        probit.like <- function (beta) {
          eta <- X %*% beta   ## 선형함수
          p <- pnorm(eta)   ## cdf 확률
          return(sum((1 - y) * log(1 - p) + y * log(p))) ## 로그우도
        }
        
        # 베타에 대한 gradient를 따로 정의해 넣어 주면 MLE를 보다 빨리 수행 
        probit.gr <- function (beta) {
          mu <- X %*% beta   ## 선형함수
          p <- pnorm(mu)  ## 확률
          u <- dnorm(mu) * (y - p) / (p * (1 - p))   ## 체인규칙
          return(crossprod(X, u))
        }
        
        formula <- low~age+as.factor(race)+smoke
        ols <- lm(formula, data=birthwt, y=TRUE, x=TRUE)
        y <- ols$y
        X <- ols$x
        
        fit <- optim(ols$coef, probit.like, gr = probit.gr, 
        						 control = list(fnscale = -1), 
                     method = "BFGS", hessian =  TRUE)
        fit.glm <- glm(formula, data=birthwt, 
                       family = binomial(link = "probit"))
        
        # par(family = "AppleGothic")
        plot(fit$par, coef(fit.glm), xlab="optim 추정치", ylab="glm 추정치")
        abline(a=0, b=1, col="red", lty=3)

      • 왈드 검정과 우도비 검정
        • 왈드검정

          최대 우도 추정치가 영가설의 값으로부터 얼마나 극단적인지를 표준오차를 통해 아래와 같이 왈드 통계량을 표현할 수 있다. t-분포를 사용하며, 자유도는 X의 데이터의 수에 feature의 수를 뺀 값이다.

          wald statistic=θ^θ0se(θ^)wald \ statistic = {{\hat{\theta}- \theta_{0}}\over{se(\hat{\theta})}}

          	lm.like <- function(theta, y, X){ 
            n <- nrow(X)
            k <- ncol(X)
            beta <- theta[1:k] 
            sigma <- exp(theta[k + 1])
            loglike <- sum(log(dnorm(y, X %*% beta, sigma)))
            return(loglike)
          }
          
          
          set.seed(1999)
          X <- cbind(1, runif(100))
          theta.true <- c(1,-1, 1)
          y <- X%*%theta.true[1:2] + rnorm(100)
          K <- ncol(X)
          
          ## 최대우도추정
          mle <- optim(c(1,1,1), lm.like, method="BFGS", 
                       control = list(fnscale = -1), 
                       hessian=TRUE, y=y, X=X)
          
          ## 추정치 추출
          beta.mle <- mle$par[1:K]
          beta.mle.var <- solve(-mle$hessian)[1:K, 1:K]
          psi.mu <- mle$par[K+1]
          psi.var <- sqrt(solve(-mle$hessian)[K + 1, K + 1])
          
          sigma.mle <- exp(psi.mu + 0.5*sqrt(psi.var)) 
          sigma.mle.var <- sigma.mle^2*exp(sqrt(psi.var) - 1)
          coef.se <- c(sqrt(diag(beta.mle.var)), sqrt(sigma.mle.var))
          coef.mle <- c(beta.mle, sigma.mle)   
          
          ## 왈드 검정값 계산
          t.stat <- coef.mle/coef.se
          pval <- 2*(1-pt(abs(t.stat), nrow(X)-ncol(X))) # pt: t분포의 cdf
          results <- cbind(coef.mle, coef.se, t.stat, pval)
          colnames(results) <- c("coef","se","t-stat","p-value")
          rownames(results) <- c("Constant","slope","sigma2")
          print(results,digits=3)

          선형회귀분석의 MLE 예시에 대해 왈드 검정 결과는 위와 같다. MLE를 통한 파라미터(기울기와 계수)가 영가설인 0과 통계적으로 유의미한 차이를 보인다고 볼 수 있다.

        • 우도비 검정

          왈드 검정은 파라미터 개별 추정치의 통계적 유의성을 검정한다면, 우도비 검정은 F-검정과 유사하게 모형 전체의 유의성 검정에 사용하며 다음과 같이 정의된다.

          λ(y)=f(y;θ,가설1)f(y;θ,가설2)\lambda(y) = {{f(y;\theta, 가설1)}\over{f(y;\theta, 가설2)}}

          우도비는 y에 대해 가설 1의 자료 생성과정에서 발생했을 가능성과 자료 2의 자료 생성과정에서 발생했을 가능성을 비교한 것이다. 가설 1의 모수는 가설 2의 선형제한으로 표현할 수 있어야 하며, 변수 선택 시에도 많이 사용한다. Wilks Theorem에 따르면 우도비에 로그를 씌우고 2를 곱하면 카이제곱 분포를 따른다. 이를 이용해, 가설 1의 가설 2에 대한 상대적 타당성을 검정할 수 있다.

    댓글

Designed by Tistory.