0. Introduction

Point process model 혹은 exponential random graph model (ergm) 등과 같은 통계적 모형은 bayesian inference가 쉽지 않다. 그 이유는 지난 [Speeding up MCMC by Delayed Acceptance and Data Subsampling] 논문에서 밝힌 바와 비슷하게, 아래와 같이 function의 normalizing constant($z(\boldsymbol \theta)$)가 parameter의 function이기 때문에 likelihood($L(\boldsymbol \theta | \mathbf x) = h(\mathbf x | \boldsymbol \theta) / z(\boldsymbol \theta)$)를 구하는 것이 computational expensive 하기 때문이다. (Posterior distribution은 다음과 같이 쓸 수 있다)

$$\pi(\boldsymbol \theta | \mathbf x) \propto p(\boldsymbol \theta) \frac{h(\mathbf x| \boldsymbol \theta)}{z(\boldsymbol \theta)} $$

위 문제를 해결하기 위한 여러 방법 중 가장 널리 쓰이는 방법은 auxiliary variable을 generate하여 likelihood의 normalizing constant 부분을 cancel out하는 double metropolis hastings (DMH) 알고리즘, 혹은 likelihood를 구하지 않고 proposed된 parameter하에서 sample을 생성하여 이것을 accept or reject하는 likelihood-free 알고리즘인 approximate bayesian computation(ABC) 알고리즘 등이 있다.

하지만 여전히 두 알고리즘 역시 theoretical한 측면에서 guarantee 되거나 computational expensive하다는 장점이 있다. 따라서 이 논문에서 저자들은 likelihood function의 normalizing constant를 Gaussian process approximation을 통하여 approximate하는 새로운 알고리즘을 제시한다.

1. Function Emulation

논문에서 중요한 개념 중 하나인 function emulation이란 함수를 계산하는 것이 어려울 때 approximate한 함수를 사용하는 방법이다. 예를 들어 위의 likelihood function의 normalizing constant $z(\boldsymbol \theta)$를 계산하는 것이 어렵다면 이 부분을 approximate하여 likelihood 계산에 사용하는 것이다.

이 방법은 point process model 혹은 ergm 과 같이 likelihood 계산이 복잡한 문제에 대해서 기존의 DMH 등의 알고리즘보다 bayesian inference를 computationally effcient하게 할 수 있다는 것이 장점이다. 저자들이 제시한 알고리즘은 다음과 같이 크게 두 단계로 구성된다.

(1) normalizing constant을 importance sampling을 사용하여 approximate

(2) Gaussian process를 사용하여 관측되지 않은 parameter의 importance sampling estimate

2. Algorithm

p-dimensional parameter $\boldsymbol \theta \in \Theta$에 대하여, $\Theta$안의 $d$개의 particle, $\psi = (\boldsymbol \theta^{(1)}, \dots,\boldsymbol \theta^{(d)})$, 을 고려해보자. 그리고 $\tilde{\boldsymbol \theta}$를 MPLE or sample mean 등의 mle의 approximation 값이라고 하자.

그렇다면 우리는 $(Z(\boldsymbol\theta^{(1)})/Z(\tilde{\boldsymbol\theta}),\dots,Z(\boldsymbol\theta^{(d)})/Z(\tilde{\boldsymbol \theta}))$의 unbiased Monte Carlo estimate를 importance sampling을 사용하여 아래와 같이 log scale로 estimate 할 수 있다 (이것의 estimate가 필요한 이유는 MCMC algorithm으로 acceptance probability를 계산할 때 해당 부분을 계산해야 하는데, 이것이 어렵기 때문에 위의 MC estimate를 사용하는 것이다).

$$ \log \widehat Z_{IS}(\boldsymbol\theta^{(i)}) = \log\Big( \frac{1}{N} \sum_{l=1}^N \frac{h(\mathbf x_l | \boldsymbol \theta^{(i)})}{h(\mathbf x_l | \tilde{\boldsymbol \theta})}\Big)$$

각 $\mathbf x_l$은 l-th Markov chain의 마지막 샘플을 의미한다. 위의 MC estimate를 각 d개의 particle에 대하여 모두 모은 것을 $\log \hat{\mathbf Z}_{IS} = (\log \hat{Z}_{IS}(\boldsymbol\theta^{(1)}),\dots,\log \hat{Z}_{IS}(\boldsymbol\theta^{(d)}))^\prime \in \mathbb R^d$ 라고 하자. 우리는 아래와 같이 Gaussian process model을 정의할 수 있다.  ($\boldsymbol \mu$ is the mean and $\mathbf u$ is a second-order stationary Gaussian process)

$$\log\hat{\mathbf Z}_{IS} = \boldsymbol \mu + \mathbf u$$

또한 새로운 parameter $\boldsymbol \theta^* \in \Theta$에 대하여, 우리는 Gaussian process의 definition을 사용하여 아래와 같은 multivariate gaussian normal model을 정의할 수 있다. (simple linear mean trend $\boldsymbol \mu = \boldsymbol \psi \boldsymbol \beta$를 가정, $\mathbf C$ and $\mathbf c$ are the covariance functions in equation (4) in the paper, and $\tau^2 \ \& \ \sigma^2$ are the parameters of covariance function)

$$\begin{align*}\begin{bmatrix} \log \hat{\mathbf Z}_{IS}\\ \log \hat{Z}_{IS}(\boldsymbol \theta^*)\end{bmatrix} &\sim \text{MVN}\begin{pmatrix} \begin{bmatrix} \boldsymbol \psi \boldsymbol\beta \\ \boldsymbol \theta^* \boldsymbol \beta \end{bmatrix}\! ,& \begin{bmatrix} \mathbf C & \mathbf c \\ \mathbf c^\prime & \sigma^2 + \tau^2 \end{bmatrix}\end{pmatrix} \end{align*}$$

따라서 우리는 새로운 $\theta^*$의 likelihood를 매번 계산하지 않고, 주어진 $d$개의 particle set을 사용하여 만든 위 Gaussian process approximation을 통하여 $\log \hat Z(\boldsymbol \theta^*) | \log \hat{\mathbf Z}_{IS}$를 구할 수 있고, 이것을 사용하여 best linear unbiased predictor(BLUP)를 구할 수 있다. (논문 equation (6) ~ (8) 참고)

논문의 알고리즘 1을 살펴보면, Part 1에서는 위 과정 중 importance sampling approximation 부분을 수행하며, Part 2에서는 Part 1에서 만든 importance sampling approximation을 사용해 새롭게 propose되는 theta에 대한 Gaussian process approximation을 수행하여 likelihood를 우회적으로 계산 후 acceptance probability를 계산해 sample의 accept 여부를 결정한다.

또한 likelihood의 normalizing constant part가 intractable한 상황이 아니라 likelihood의 계산 자체가 intractable한 경우 역시 likelihood function 자체를 emulation 하는 방법을 통해 논문의 알고리즘 2와 같이 극복할 수 있다.

다른 중요한 한가지는 d 개의 parameter에 대한 particle을 만들 때, 각 particle이 parameter space의 중요한 부분을 잘 cover하도록 만드는 것이다. 저자들은 DMH or ABC algorithm을 사용하면 효과적으로 중요한 parameter space를 cover하는 particle을 만들 수 있음을 알고리즘 3, 4를 통하여 밝힌다.

3. Application

마지막으로 저자들이 만든 알고리즘을 기존의 DMH와 비교한 real data application 실험을 진행하였다. 실험은 총 두 가지로 첫번째는  실제 social network data와 simulated data를 ergm을 사용하여 fitting한 모델의 parameter를 estimate한 결과와 computation time을 비교한 실험이며, 두번째는 attraction-repulsion point process 모델의 Parameter 추정 및 computation 시간을 비교한 것이다.

중요한 비교 결과는 두 실험 모두에서 DMH 알고리즘에 비해 크게 computation time을 단축시켰다는 것이다. 아래에 첨부한 Figure 1을 보면 큰 데이터 셋에 대하여 computing 시간이 DMH에 비교해 훨씬 줄어든 것을 확인할 수 있으며, parameter estimation 결과는 비슷한 것을 확인할 수 있다.

Reference:

Park, J., & Haran, M. (2020). A function emulation approach for doubly intractable distributions. Journal of Computational and Graphical Statistics, 29(1), 66-77.

Link:

https://www.tandfonline.com/doi/full/10.1080/10618600.2019.1629941?casa_token=XKpfYl743TUAAAAA%3AfFdKaRxkgDgl68nTbjn0QvRtI1ExBBrtobCw_eb707hZwVwM1Pr1m2cCAJMg7X0OVimOndrthvTExw

0. Introduction

Bayesian framework는 MCMC나 ABC 등과 같은 simulation 기반 실험에서 사용된다. 이는 likelihood function을 구하지 않고 posterior 분포를 구할 수 있다는 장점을 가지고 있으나, 고차원의 세팅에서는 수렴 속도가 느리거나 posterior sampling 결과가 좋지 않다는 단점을 가지고 있기도 하다.

이 논문의 저자들은 이러한 점에서 문제를 제기하고, complex (intractable) model에서 neural network를 사용하여 관심있는 posterior 분포의 파라미터들을 구하는 방법을 제기한다. 이 방법의 장점은 posterior를 target으로 설정하여 computational cost를 줄일 수 있고 고차원의 데이터를 neural network로 다루기 용이하다는 것 등이 있다.

1. Model - VaNBayes

이들이 사용한 bayesian nerual network의 이름은 VaNBayes이며 몇가지 notation을 살펴보자. 첫째로 $\mathbf{Y} = (Y_1,\dots,Y_n)^\top$ 은 $f(\mathbf{Y} \mid \mathbf{\theta})$를 likelihood로 가지는 response vector 이다. Model parameter는 $\mathbf{\theta} = (\theta_1,\dots,\theta_p)^\top$라고 하자. 이때 Bayesian model을 다음과 같이 쓸 수 있다.

\begin{equation} \mathbf{Y}\mid\mathbf{\theta} \sim f(\mathbf{Y} \mid \mathbf{\theta}) \end{equation}

\begin{equation}\mathbf{\theta} \sim \pi(\mathbf{\theta}) \end{equation}

$\pi$는 $\boldsymbol{\theta}$의 prior 분포이며 $f$ 또는 $\pi$가 intractable함을 가정한다. 또한 parameter $\boldsymbol{\theta}$에 대한 Q-dimensional summary의 추론에 관심있음을 가정하며 이를 $\boldsymbol{\gamma} = (\gamma_1,\dots,\gamma_Q) = G(\boldsymbol{\theta})$로 표기한다.

그 후 우리가 관심있는 $\boldsymbol{\gamma}$의 posterior는 다음과 같이 정의한다.

$$ \boldsymbol{\gamma}\mid\mathbf{Y} \sim p\{\boldsymbol{\gamma}\mid a(\mathbf{Y}; \mathbf{W})\}$$

분포 $p$는 사전에 알려져있거나 실험자에 의해 설정되며, $a(\mathbf{Y} ; \mathbf{W}) = \{a_1(\mathbf{Y} ; \mathbf{W}_1),\dots,a_J(\mathbf{Y};\mathbf{W}_J)\}, \mathbf{W} = \{\mathbf{W}_1,\dots,\mathbf{W}_J\}$에서 각 $a_j$는 weight를 $\mathbf{W}_j$로 가지는 neural network이다. 각 weight는 $\boldsymbol{\gamma}$에 대한 marginal posterior distribution 정보를 담도록 학습되어야 한다. 또한 $\mathbf{Y}$를 그대로 사용하기 보다는 차원을 축소하기 위해 summary statistic이 존재하는 경우 summary statistic $\mathbf{Z}=S(\mathbf{Y})$를 대신 사용하기도 한다.

요약하자면 parameter를 그대로 inference하는 것이 목표가 아닌, parameter에 대한 Q-dimensional summary의 inference를 목표로 하고 있으며 (차원 축소의 효과), 그 수단으로 neural network를 사용한다는 것을 알 수 있다.

VaNBayes를 실행하기 위해 필요한 다른 것은 training distribution $\Pi$이며 $\boldsymbol{\theta}$를 generate하는 역할을 하기 때문에 prior distribution인 $\pi$와 같은 support를 가져야 한다. 알고리즘은 training distribution에서 $\boldsymbol{\theta}$를 sampling하는 것에서부터 시작하며 다음과 같이 쓸 수 있다.

$$\boldsymbol{\theta}_i  \overset {iid} \sim\Pi(\boldsymbol{\theta}) \quad \text{and} \quad \boldsymbol{Y}_i \mid \boldsymbol{\theta}_i \overset {indep} \sim f(\boldsymbol{Y}\mid\boldsymbol{\theta}_i) \ \ \text{for } \ i \in \{1,\dots,N\}$$

각 posterior distribution에 곱해지는 weight인 $w$는 $w_i = \pi(\boldsymbol{\theta}_i) / \Pi(\boldsymbol{\theta}_i)$로 쓸 수 있다. 여기까지 VaNBayes를 위한 사전 준비는 모두 끝났다. VaNBayes 알고리즘은 다음과 같다.

앞서 설명한 training distribution, likelihood 등이 필요하며, 우선 training distribution에서 parameter를 sampling하며 이것을 바탕으로 새로운 data $Y_i$를 generate 하며, 이것의 요약 통계량 $Z$와 parameter theta의 요약 값인 gamma를 계산한다. (1~4)

그 후 neural network의 weight인 $\boldsymbol{W}$를 조정해가며 log-posterior distribution의 likelihood가 최대가 되는 neural network weight $\widehat{\mathbf{W}}$를 구한다. (5)

최종적으로 true posterior에 근사하게 되는 neural network의 weight기반 estimated posterior를 구할 수 있다. (6) 이때 우리는 theta에 대한 posterior가 아닌 gamma에 대한 posterior를 구하는 것에 관심이 있으므로 gamma에 대한 posterior를 도출해낸다.

VaNBayes가 작동하는 이론적 이유는 논문 2.3에 잘 나와있다. true joint posterior와 marginal * neural network posterior 사이의 Kullback-Leibler divergence가 최소화되는 방향으로 neural network weight가 학습된다는 것이 주된 내용이다.

또한 실험을 통해서도 도출된 posterior가 실제에 잘 근사하는지 평가하는 과정이 있다. log score $LS = \sum_{v=1}^V \log p\{\boldsymbol{\gamma}_v \mid a(\boldsymbol{Z}_v ; \widehat{\mathbf{W}}) \}$ 를 사용하여 서로다른 data $(\boldsymbol{\gamma}_v, \boldsymbol{Z}_v)$에 대한 특정 model $a(\cdot)$의 fit을 평가한다.

2. 실험

마지막으로 VaNBayes를 사용한 다양한 실험과 결과 등에 대해 살펴보자.

2.1 Sparse linear regression

- 실험 세팅 및 순서

먼저 아래와 같은 Sparse linear regression model에서 coefficient beta를 잘 capture하는지 실험하고, 최종적으로 그 결과를 MCMC를 사용해 sampling한 것과 비교하였다.

최종 목표는 $\sigma$와 posterior inclusion probability $\text{PIP}_j = \text{Prob}(\beta_j \neq 0 \mid \mathbf{Y}) \ \text{for } j \in \{1,\dots,p\}$를 estimate하는 것이다.

$$Y_i = \beta_0 + \sum_{j=1}^p X_{ij} \beta_j + \epsilon_i$$

$\epsilon_i \overset{iid} \sim N(0, \sigma^2), i \in \{1,\dots,n\}$, 각 X들은 평균 0, 분산 1이며 $Cor(X_{ij},X_{ik})=\rho^{\lvert j-k \rvert}$인 상관관계를 가진다.

우리가 관심있는 model parameter인 $\beta$는 true 값이 $\beta_0=0, \beta_1=\beta_2=\beta_6=0.5$이며 1, 2, 6을 제외한 나머지 $\beta_j=0$ 이다. $\rho=0.5, \sigma=1$로 설정했으며 각 50개의 데이터에 대한 총 100번의 실험을 반복하였다. covariate $p$의 개수는 10개, 20개인 경우를 각각 진행하였다.

전반적인 실험의 순서는 아래와 같다.

1. 위 모델에서 data generation (100 sets)

2. summary statistic과 gamma 사용하여 neural network training

3. PIP와 sigma estimate

앞서 살펴본 VaNBayes를 위해 필요한 것이 $\gamma, Z$ 임을 알았다.

우선 $\gamma$의 경우 $\gamma_j = \mathbb{I}(\beta_j \neq 0) \in \{0,1\}, j \in \{1,\dots,p\}, \gamma_{p+1} = \sigma, \gamma_{p+1+i} =Y_{n+i}, i\in \{1,\dots,10\}$ 로 설정하였다.

그 다음 $\mathbf{Z}$를 설정한 방법에 대해 살펴보면, $(\beta_0,\dots,\beta_p)$의 least squares estimate & residual standard deviation & standard deviation of the least square estimates의 총 $p+3$의 summary statistic을 사용하였다.

또한 두 개의 hidden layer를 가지는 feed forward neural network $a_j(\mathbf{Z} ; \mathbf{W}_j)$ 를 구성하여 logistic regression model $logit\{P(\gamma_j=1 \mid \mathbf{Z})\} = a_j(\mathbf{Z} ; \mathbf{W}_j)$를 fit 하였다.

- 실험 결과

위에서 적은 내용 이외에도 neural network weight의 node수를 정하는 방법 등이 첨부되어 있으며, PIP 를 예측한 실험 결과 및 MCMC 결과와의 비교는 아래 그림과 같다.

covariate의 수가 증가하는 경우에 VaNBayes가 일부 PIP를 underestimate하는 경향이 있으나 전반적으로는 MCMC와 비슷한 예측을 하는 것을 알 수 있다. 또한 아래 그림처럼 Sigma 역시 MCMC의 결과와 비슷하게 예측한 것을 알 수 있다.

2.2 Autologistic regression model

다음으로 살펴볼 실험은 autologistic regression model에서의 parameter beta의 posterior의 estimate 실험이다. 간단히 말해서 autologistic은 covariate $X_i$와 node i 의 이웃들을 고려하여 $Y_i\in \{0,1\}$의 값을 가지게 된다. 주변 데이터와의 dependency를 모델링 한 것이라고 볼 수 있겠다.

autologistic 실험에서도 앞선 실험에서 했던 것과 비슷한 세팅이 필요하다.

- 실험 세팅 및 순서

autologistic model의 full conditional distribution of one node given the others는 다음과 같다.

$$logit\{\text{Prob}(Y_i=1\mid Y_j, j \neq i)\} = logit(\kappa_i) + \phi \sum_{j\in \mathcal{N}_i} (Y_j - \kappa_j)$$

$\mathcal{N}_i$는 node i와 연결된 다른 노느들의 indices의 집합이다. 또한 $logit(\kappa_i) = \mathbf{X}_i^\top \boldsymbol{\beta}$는 covariate effect $\boldsymbol{\beta} = (\beta_1,\dots,\beta_p)$를 가지는 logistic regression probability이며 $\phi>0$은 dependence의 강도를 나타낸다.

관심있는 parameter 역시 $\boldsymbol{\beta}, \log(\phi)$이며 이들의 prior는 standard normal distribution이다. neural network를 train하기 위한 summary statistic $\mathbf{Z}$로 사용한 것은 각 $\beta_j$들의 logistic regression estimates & $\widehat{\boldsymbol{\beta}}_{GLM}$ & Geary's C statistics이다.

$\boldsymbol{\gamma}$로는 parameter를 그대로 사용하여 $\boldsymbol{\gamma}=\boldsymbol{\theta}=[\beta_1,\dots,\beta_p,\log(\phi)]$ 이다.

20 by 20 grid에서 400개의 data를 simulation하였고, covariate effect의 수(p) = 5로 설정하였다. 이것을 100,000번 반복하여 neural network를 train하였다. 그 결과는 아래와 같은데 dependency,$\phi$,가 작은 경우는 posterior median이 조금 벗어나있긴 하지만 다른 parameter들에 대해서는 true parameter인 red dot과 비교했을때 매우 훌륭한 예측을 한 것을 볼 수 있다.

3. 요약

지금까지 bayesian neural network를 사용하여 parameter를 estimation하는 방법에 대해 살펴보았다. 여기서 주된 내용은 parameter에 대한 direct한 estimation이 아닌 최종적으로 posterior distribution을 estimate하는 neural network를 만들었다는 점이다.

또한 parameter와 그에 기반해 생성한 data $Y$를 그대로 input하여 사용하기 보다는 차원을 축소하여 computational efficiency를 높이기 위해 summary statistic과 같은 input을 사용하였으며 Q-dimensional summary $\boldsymbol{\gamma}$에 대한 inference를 진행하였다.

알고리즘을 간략히 요약하자면 최종 목표가 observed data $\mathbf{Y}_0$의 parameter에 대한 estimation이라고 했을때, VaNBayes에서는 training distribution에서 parameter를 sampling하고, 그를 기반해서 새로운 $\mathbf{Y}_j$를 sampling한다.

그 후에 parameter와 data의 차원을 줄인 후 posterior distribution의 log-likelihood를 최대로 만드는 neural network weight $\widehat{\mathbf{W}}$를 찾는다.

Ref: Maceda, E., Hector, E. C., Lenzi, A., & Reich, B. J. (2024). A variational neural Bayes framework for inference on intractable posterior distributions. arXiv preprint arXiv:2404.10899.