비대칭 제안서 배포를 통한 대도시 해 스팅 이해 있습니다. 생성 Yt∼q(y|xt)Yt∼q(y|xt)Y_t \sim q(y|x^t) Xt+1={Yt,xt,with probabilityρ(xt,Yt),with

모델의 매개 변수를 추정하기위한 코드를 작성하기 위해 Metropolis-Hastings 알고리즘을 이해하려고 노력했습니다 (예 :

f(x)=a∗x

). 참고 문헌에 따르면 Metropolis-Hastings 알고리즘에는 다음 단계가 있습니다.

  • 생성
    Yt∼q(y|xt)


  • Xt+1={Yt,with probabilityρ(xt,Yt),xt,with probability1−ρ(xt,Yt),

ρ(x,y)=min(f(y)f(x)∗q(x|y)q(y|x),1)

몇 가지 질문을하는 방법 :

  • 참고 문헌은 가 대칭 분포이면 비율 가 1이되고 알고리즘을 메트로폴리스라고합니다. 그 맞습니까? Metropolis와 Metropolis-Hastings의 유일한 차이점은 첫 번째가 대칭 분포를 사용한다는 것입니다. “Random Walk”Metropolis (-Hastings)는 어떻습니까? 다른 두 가지와 어떻게 다릅니 까?q ( x | y ) / q ( y | x )
    q

    q(x|y)/q(y|x)

  • 온라인에서 찾은 대부분의 예제 코드는 가우스 제안 분포 를 사용하므로 여기서 는 가능성 비율입니다. 제안 배포가 포아송 배포 인 경우 어떻게합니까? 비대칭 분포를 사용할 때 그 비율이 1이 아닌 이유를 합리적으로 이해한다고 생각하지만 수학적으로 이해하거나 코드로 구현하는 방법을 잘 모르겠습니다. 누군가가 비대칭 제안 분포를 사용하는 Metropolis-Hastings 알고리즘의 간단한 코드 (C, python, R, pseudo-code 또는 원하는 코드) 예제를 제공 할 수 있습니까?ρ ( x , y ) = 최소 ( f ( y ) / f ( x ) , 1 ) f ( y ) / f ( x )
    q

    ρ(x,y)=min(f(y)/f(x),1)

    f(y)/f(x)

답변

참고 문헌은 q가 대칭 분포이면 비율 q (x | y) / q (y | x)가 1이되고 알고리즘을 메트로폴리스라고합니다. 그 맞습니까?

예, 맞습니다. 메트로폴리스 알고리즘은 MH 알고리즘의 특수한 경우입니다.

“Random Walk”Metropolis (-Hastings)는 어떻습니까? 다른 두 가지와 어떻게 다릅니 까?

무작위 보행에서, 제안 분배는 체인에 의해 마지막으로 생성 된 값에서 각 단계 후에 중심이 조정됩니다. 일반적으로 랜덤 워크에서 제안 분포는 가우시안이며이 경우이 랜덤 워크는 대칭 요구 사항을 충족하며 알고리즘은 대도시입니다. 비대칭 분포로 “의사”랜덤 보행을 수행 할 수 있다고 가정하면 제안이 기울기의 반대 방향으로 너무 표류하게됩니다 (왼쪽으로 치우친 분포는 오른쪽으로 제안을 선호합니다). 왜 당신이 이것을 할 것인지 잘 모르겠지만, 당신은 할 수 있고 대도시 헤이스팅스 알고리즘 일 것입니다 (즉, 추가 비율 항이 필요합니다).

다른 두 가지와 어떻게 다릅니 까?

비임의 보행 알고리즘에서 제안 분포는 고정되어 있습니다. 랜덤 워크 변형에서는 제안 반복의 중심이 각 반복마다 변경됩니다.

제안 배포가 포아송 배포 인 경우 어떻게합니까?

그런 다음 대도시 대신 MH를 사용해야합니다. 아마도 이것은 불연속 분포를 샘플링하는 것일 수 있습니다. 그렇지 않으면 불연속 함수를 사용하여 제안서를 생성하지 않을 것입니다.

어쨌든 샘플링 분포가 잘 리거나 기울어 짐에 대한 사전 지식이있는 경우 비대칭 샘플링 분포를 사용하려고하므로 대도시 헤이스팅스를 사용해야합니다.

누군가 나에게 간단한 코드 (C, python, R, pseudo-code 또는 원하는 코드)를 줄 수 있습니까?

대도시는 다음과 같습니다.

Metropolis <- function(F_sample # distribution we want to sample
                      , F_prop  # proposal distribution
                      , I=1e5   # iterations
               ){
  y = rep(NA,T)
  y[1] = 0    # starting location for random walk
  accepted = c(1)

  for(t in 2:I)    {
    #y.prop <- rnorm(1, y[t-1], sqrt(sigma2) ) # random walk proposal
    y.prop <- F_prop(y[t-1]) # implementation assumes a random walk.
                             # discard this input for a fixed proposal distribution

    # We work with the log-likelihoods for numeric stability.
    logR = sum(log(F_sample(y.prop))) -
           sum(log(F_sample(y[t-1])))

    R = exp(logR)

    u <- runif(1)        ## uniform variable to determine acceptance
    if(u < R){           ## accept the new value
      y[t] = y.prop
      accepted = c(accepted,1)
    }
    else{
      y[t] = y[t-1]      ## reject the new value
      accepted = c(accepted,0)
    }
  }
  return(list(y, accepted))
}

바이 모달 분포를 샘플링하기 위해 이것을 사용해 봅시다. 먼저, 우리의 propsal에 임의의 보행을 사용하면 어떻게되는지 봅시다 :

set.seed(100)

test = function(x){dnorm(x,-5,1)+dnorm(x,7,3)}

# random walk
response1 <- Metropolis(F_sample = test
                       ,F_prop = function(x){rnorm(1, x, sqrt(0.5) )}
                      ,I=1e5
                       )
y_trace1 = response1[[1]]; accpt_1 = response1[[2]]
mean(accpt_1) # acceptance rate without considering burn-in
# 0.85585   not bad

# looks about how we'd expect
plot(density(y_trace1))
abline(v=-5);abline(v=7) # Highlight the approximate modes of the true distribution

여기에 이미지 설명을 입력하십시오

이제 고정 제안서 배포를 사용하여 샘플링을 시도하고 어떻게되는지 살펴 보겠습니다.

response2 <- Metropolis(F_sample = test
                            ,F_prop = function(x){rnorm(1, -5, sqrt(0.5) )}
                            ,I=1e5
                       )

y_trace2 = response2[[1]]; accpt_2 = response2[[2]]
mean(accpt_2) # .871, not bad

처음에는 괜찮아 보이지만, 만약 우리가 후부 밀도를 살펴보면 …

plot(density(y_trace2))

여기에 이미지 설명을 입력하십시오

로컬 최대 값에 완전히 고정되어 있음을 알 수 있습니다. 우리가 실제로 제안서 배포를 중심으로 했으므로 이것은 놀라운 일이 아닙니다. 우리가 이것을 다른 모드에 집중 시키면 같은 일이 일어납니다 :

response2b <- Metropolis(F_sample = test
                        ,F_prop = function(x){rnorm(1, 7, sqrt(10) )}
                        ,I=1e5
)

plot(density(response2b[[1]]))

두 가지 모드 사이 에서 제안을 제거 할 수 있지만 두 가지 모드 중 하나를 탐색 할 수 있도록 분산을 실제로 높게 설정해야합니다.

response3 <- Metropolis(F_sample = test
                        ,F_prop = function(x){rnorm(1, -2, sqrt(10) )}
                        ,I=1e5
)
y_trace3 = response3[[1]]; accpt_3 = response3[[2]]
mean(accpt_3) # .3958!

제안서 배포의 중심 선택이 어떻게 샘플러의 수락 률에 큰 영향을 미치는지 살펴보십시오.

plot(density(y_trace3))

여기에 이미지 설명을 입력하십시오

plot(y_trace3) # we really need to set the variance pretty high to catch
               # the mode at +7. We're still just barely exploring it

우리는 여전히 두 가지 모드에 가까워지고 있습니다. 이것을 두 모드 사이에 직접 놓아 봅시다.

response4 <- Metropolis(F_sample = test
                        ,F_prop = function(x){rnorm(1, 1, sqrt(10) )}
                        ,I=1e5
)
y_trace4 = response4[[1]]; accpt_4 = response4[[2]]

plot(density(y_trace1))
lines(density(y_trace4), col='red')

여기에 이미지 설명을 입력하십시오

마지막으로, 우리는 찾고 있던 것에 더 가까워지고 있습니다. 이론적으로, 만약 우리가 샘플러를 충분히 오래 운영하게한다면, 우리는이 제안서 분포에서 대표 샘플을 얻을 수 있지만, 랜덤 워크는 유용한 샘플을 매우 빨리 생성했으며, 어떻게 후자가 추정되었는지에 대한 지식을 이용해야했습니다. 고정 된 샘플링 분포를 조정하여 유용한 결과를 얻습니다 (실제로 아직 알려지지 않은 결과 y_trace4).

나중에 대도시 헤이스팅의 예로 업데이트하려고합니다. 메트로폴리스 헤이스팅스 알고리즘을 생성하기 위해 위의 코드를 수정하는 방법을 상당히 쉽게 볼 수 있어야합니다 (힌트 : 추가 비율을 logR계산 에 추가하면 됩니다).

답변