이 글에서 PRML의 11장에서 언급한 마르코프 체인 몬테카로(이하 MCMC)는 구체적으로python에서 Metropolis algorithm를 실현하는 것이다.
대응하는 주피터 노트북은 필자github 저장소에 있다.전체적인 연재 방침과 PRML의 다른 알고리즘의 실현에 대해서는 보십시오연재의 총결산 페이지.
지난번 10장 보도와는 달리 이번 이론과 설치는 모두 간단한 보도일였다.

1 이론 복습

1.1 설정

주어진 집합 이 $X$의 확률 분포 (확률 밀도 함수) $p$를 가정합니다.
이후의 편의를 위해 확률 분포(확률 밀도 함수) $p$는 $p(x) =\tilde{p}(x)/Z_와 같다p$를 쓸 수 있다고 가정합니다.여기, $\tilde{p}(x)$는 이미 알고 있지만, 표준화 상수 $Z_만약 p$가 알 수 없다고 가정해도 삼 가능합니다.
MCMC(PRML11장에서 기술한 샘플링 알고리즘에만 국한되지 않음)의 목표는 $p$에 따라 $X$점의 샘플링을 생성하는 것이다.
MCMC는 마르코프 체인을 사용하여 이 샘플을 생성합니다.

1.2 Metropolis-Hastings algorithm

이 글에서 저는 먼저 MCMC의 방법 중 하나인 Metropolis-Hastings algorithm를 소개하고 싶습니다.
우선, 제안 분포 $q$를 결정합니다.$q(x｜x')$는 현재 상태가 $x'$일 때 $x$를 다음 상태의 후보자로 받을 확률을 나타냅니다. $만약 q(\cdot\ux')$가 쉽게 샘플링된다고 가정합니다.또한 $x, x'$가 비대칭이면 확률 (확률 밀도) $q(x｜x') $의 값을 쉽게 평가할 수 있다고 가정합니다.
이 전제에서 Metropolis-Hastings algorithm은 다음과 같은 작업을 수행합니다.

현재 상태의 $x^ {current}$(첫 번째 단계라면 외부에서 제시한 초기 상태의 $x^ {ini}$), $q(\_cdot_x^ {current}) $에 따라 다음 상태의 후보 $x^ {tmp}$사를 생성합니다.

아래 상태 확인 $x^ {next}$
$$
\begin{align}
x^{next} =
\begin{cases}
x^{tmp} & (\mbox{with probability } A(x^{tmp}, x^{current}))\\
x^{current} & (\mbox{with probability } 1 - A(x^{tmp}, x^{current}))
\end{cases},
\end{align}
$$
하지만
$$
\begin{align}
A(x,y) :=\min\left( 1,\\frac{\tilde{p}(x)}{\tilde{p}(y)}\cdot\frac{q(y|x)}{q(x|y)}\right)
\end{align}
$$
현재 상태가 $y$일 때 다음 상태의 후보 $x$를 받아들일 확률오을 나타냅니다.

이렇게 마르코프 체인을 정의하면 $p$가 변하지 않고 분포됩니다.이것은 상세한 균형이 성립된 것을 볼 수 있다육.
(비슷하게) 독립된 샘플을 원한다면, 이렇게 얻은 샘플 사이의 간격 칠: $x^ {(0)}, x^ {(1)},\dots$를 기록해야 합니다.
이론편은 여기서 끝내자!

두 공식에서 코드로

그럼 공식을 코드로 표시합시다.그래도 이번에는 담백하다.
이 절은 MetropolisSampler류를 실현할 것이다.
다음은 $X$가 $D$비오클리드 공간이라고 가정하고, 제안된 분포는 협방차 $\Sigma$의 $D$비고스 분포입니다.
$$
\begin{align}
q(x|y) =\frac{1}{(2\pi)^{D/2}\sqrt{\det\Sigma}}\exp\left[ -\frac{1}{2}(x-y)^T\Sigma^{-1} (x-y)\right]
\end{align}
$$
사용자 정의 모양새를 정의합니다.
또한 $q(x｜y)$는 $x, y$의 교체가 대칭적이기 때문에 이 경우 수리 확률 $A$는 더욱 간결하게 쓸 수 있다.
$$
\begin{align}
A(x,y) :=\min\left( 1,\\frac{\tilde{p}(x)}{\tilde{p}(y)}\right)
\end{align}
$$

2.1 데이터 속성 및 방법

MetropolisSampler 클래스에는 다음과 같은 데이터 속성과 방법이 있습니다.

2.1.1 데이터 속성

D : $D$. i.e., 사고공간 $X$의 차원.

sigma: $\Sigma$, i.e., 분포된 협방차 행렬을 권장합니다.

p: $\tilde{p}$, 샘플을 찾으려는 확률(일정 표준화되지 않음)의 확률 밀도 함수를 나타냅니다.

2.1.2 방법

_A: $A$, 수리 확률을 계산하는 함수입니다.

sample: Metropolis algorithm에 따라 샘플을 채취하여 샘플의 함수를 되돌려줍니다."각 레코드 샘플"을 지정할 수 있습니다.

2.2 코드

class MetropolisSampler:
    def __init__(self, D, sigma, p):
        self.D = D # the dimension of the output space       
        self.sigma = sigma # the covariance matrix of the proporsal kernel
        self.p = p # the density function of the target distribution

    def _A(self, x, y):
        '''
        The method which returns the acceptance probability, 
        depending on the current state y and the candidate state x

        Parameters
        ----------
        x : 1D numpy array
            1D numpy array representing the candidate state
        y : 1D numpy array
            1D numpy array representing the current state
        Returns 
        ----------
        A : float
            The acceptance probability A(x, y), 
            where y is the current state, and x is the candidate state.
        '''
        denom = self.p(y)  # for avoiding zero division error
        if denom == 0:
            return 1.0
        else:
            return min( 1.0, self.p(x) / denom )

    def sample(self, xini, N, stride):
        '''
        The method that performs sampling using Metropolis algorithm, and returns the samples

        Parameters
        ----------
        xini : 1D numpy array
            1D numpy array representing the initial point
        N : int
            An integer representing the number of samples required.
        stride : int
            An integer representing how frequent samples are recorded.
            More specifically, we only record sample once in `stride` steps

        Returns 
        ----------
        X : 2D array
            (N, self.D) array representing the obtained samples, where X[n] is the n-th sample.
        '''
        X = np.zeros((N, self.D))  # array for recording the sample
        x = xini
        cnt = 0
        for i in range((N-1) * stride + 1):
            xtmp = np.random.multivariate_normal(x, self.sigma)  # sample from the proporsal kernel (gaussian )
            tmprand = random.random()
            if tmprand < self._A(xtmp, x):
                x = xtmp
            if i % stride == 0:
                X[cnt] = x
                cnt += 1
        return X

3 실험

그러면 이MetropolisSampler로 실제 샘플을 생성합시다.얻은 직사각형도와 평균/방차의 값과 확률 밀도의 도표와 평균/방차를 비교한다.
여기서 구체적으로 다음과 같은 두 가지 분포를 고려한다.

가마 분포(1차원)

2차원 정적 분포

3.1γ분포

$$
\begin{align}
& p(x) =\frac{1}{\Gamma(k)\theta^k} x^{k-1} e^{-\frac{x}{\theta}}\\
&\mbox{mean} : k\theta\\
&\mbox{variance} : k\theta^2
\end{align}
$$
실험에서, 우리는 $k=3,\theta=1$를 고려할 것이다.sigma=np.array([[1.0]]), xini=np.array([1.0]), N=5000, stride=10:

sample mean : 3.0866187342766973
sample variance : 3.284016881941356

3.2 2차원 고스 분포

$$
\begin{align}
& p(x) =\frac{1}{2\pi\sqrt{\det\Sigma}}\exp\left( -\frac{1}{2} x^T\Sigma^{-1} x\right) \\
&\Sigma^{-1} =\begin{pmatrix}
10 & -6 \\
-6 & 10
\end{pmatrix}\\
&\Sigma =
\frac{1}{32}
\begin{pmatrix}
5 & 3 \\
3 & 5
\end{pmatrix}
=\begin{pmatrix}
0.15625 & 0.09375 \\
0.09375 & 0.15625
\end{pmatrix}
\end{align}
$$

sample mean : [0.0088419  0.00923985]
sample covariance : [[0.16297775 0.09611609]
 [0.09611609 0.15594885]]

4 총결산

본고는 Metropolis-Hastings algorithm를 조직하고 Metropolis algorithm를 바탕으로 하는 샘플링을 실현했다.
실험 결과에 의하면 견본의 직사각형도와 확률 밀도는 눈이 보는 범위 내에서 큰 모순이 없고 평균/방차도 이론이 얻은 것과 큰 편차가 없다.
그러나 샘플링법에 대해 나는 알고리즘을 실시하는 것 자체가 대화의 일부분일 뿐이고 이보다 더 실천하기 어려운 부분이라고 생각한다.
갑자기 생각만 해도

평균치나 방차의 정확한 값을 모를 때 오차를 어떻게 평가해야 합니까?

샘플 수를 늘릴 때 오차는 어떻게 표현합니까?

샘플 간에 비어 있는 절차가 커질 때 샘플 간의 관련성은 어떻게 표현합니까?

더 복잡한 모델로 주행할 때 어떤 행동을 볼 수 있을까?

잠깐만.이 점에 관해서 필자는 결코 이해하지 못하는데, 이 문장에서는 닿을 수 없다.이것들에 관해서 나는 앞으로 공부할 과제로 삼고 싶다.
다만, 제가 현 상황에 대한 이해에 따라 발굴이 부족한 점도 있습니다.자세한 내용은 요약된 섹션을 참조하십시오.  ↩
정확히 말하면 측정 가능한 공간(가능성)은 적당한 $\sigma$대수를 정의한다.  ↩
전형적인 예는 매개 변수의 후분포를 고려할 수 있습니다. 매개 변수를 $\theta$로, 데이터를 $\mathcal {D}$로, 유사한 함수 $p(\mathcal {D}\theta)$와 미리 분포된 $p(\theta)$를 제시합니다.이 모델에서 데이터를 관측할 때 매개 변수의 후방 분포는 $p(\\mathcal {D})=p(\mathcal {D}\theta)p(\theta)/p(\mathcal {D})$입니다.분자 $p (\mathcal {D}\theta) p (\theta) $는 모델에서 쉽게 이해할 수 있지만, 분모 $p (\mathcal {D}) $는 많은 상황에서 계산하기 어렵다.  ↩
여기서 "$q(\cdot\x)$에서 샘플링하기 쉽다"는 가정은 유효합니다.  ↩
$A$를 계산하는 부분에서'$q(x｜x')$는 계산하기 쉽다'는 가설이 유효합니다.그러나 만약 $q$가 두 매개 변수에 대해 대칭적이라면 분자 분자로 상쇄하기 때문에 이 가설은 필요 없다.  ↩
나는 PRML(11.45)만 계산하는 것은 부족하다고 생각한다..."후보자의 샘플이 수리되지 않았을 때 앞의 상태가 이렇게 다음 상태로 바뀌었다"는 것이다.이 과정을 포함하여 이전 확률은 $T(z,z)=q(z)｜z(A(z),z)+\delta(z-z)\intdz'q(z)｜z)\left[1-A(z),z)\right]$로 표시할 수 있다.이 변환 확률과 $p$의 상세한 균형이 성립된 것을 쉽게 볼 수 있습니다.  ↩
얼마를 비우는 게 좋을지 신경 쓰이지만 앞으로의 과제로 삼고 싶어요.직관적으로 말하면 더 높은 모멘트를 계산하고 싶을 때일수록 간격을 넓힐 필요가 있지만 정량을 조사하지 않으면 어떻게 되는지 알 수 없다.  ↩