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

대부분의 MH algorithm은 단순한 모델에서는 적용이 잘 되지만, full data가 너무 크거나 intractable한 constant term을 포함하고 있어 이것의 likelihood를 얻기 어려운 경우는 적용이 어렵다는 단점이 있다. 예를 들어 Ising model과 같이 intractable constant를 포함한 모델에서 full data를 사용하여 likelihood를 계산하는데 필요한 computation이 grid 수가 증가함에 따라 2의 제곱배로 늘어난다.

따라서 이를 극복하기 위해 다양한 MH 알고리즘이 개발되어 왔다. Double Metropolis Hastings (DMH)나 Approximate Bayesian Computation (ABC)도 그것의 예시들이다. 이 논문에서 다루는 알고리즘은 Delayed Acceptance MH (DAMH) 알고리즘으로, 간단히 말하자면 두 단계로 나누어서 sample의 accept / reject 를 결정하는 알고리즘이다. 첫번째 단계에서는 random subsample을 기반으로 likelihood function을 approximate하여 잠재적인 sample의 accept 여부를 결정한다. 이 단계에서 accept된 sample에 한하여 두번째 단계로 넘어가게 되는데, 두번째 단계에서는 full data의 likelihood를 계산하여 최종적인 sample의 수락 여부를 결정한다.

이렇게 되면 sample을 propose한 후 매번 full data의 likelihood를 계산할 필요 없이, accept될 가능성이 보이는 sample에 대해서만 likelihood of full data를 계산하면 되기 때문에 기존의 MH 알고리즘보다 computation 측면에서 더욱 효율적이다.

하지만 앞서 말했듯 full data의 likelihood를 계산하는 자체가 불가능하거나 어려운 경우가 많다. 따라서 이들은 subsample을 사용한 likelihood의 estimation 방법을 제안하여 DAMH의 기존 어려움을 해결하고자 하였으며 이들의 방법이 기존의 DAMH 방법보다 더욱 효율적이고 계산 속도를 높일 수 있음을 보였다.

1. Likelihood estimator

k-th observation의 log density의 notation은 다음과 같으며 아래에서 $(y_k, x_k, \theta)$는 각각 k-th response, covariates, and model parameter를 의미하며 data size는 총 $n$개이다.

$$l_k(\theta) = \log p(y_k | \theta, x_k)$$

또한 위의 log density를 사용하여 다음과 같이 likelihood function을 정의하며, $l(\theta) = \sum_{k=1}^n l_k(\theta)$는 full data의 likelihood function을 나타낸다.

$$p(y|\theta) = \exp [l(\theta)]$$

위의 likelihood function을 estimate하기 위하여 size $m$의 sample을 사용한다. 하지만 Simple random sampling (SRS) 방법을 사용하여 sampling을 하기 보다는 이들은 control variates 방법을 활용하여 sampling을 수행한다. 먼저 다음과 같이 log likelihood를 decompose 한다.

$$\begin{align} l(\theta) & = \sum_{k=1}^n q_k(\theta) + \sum_{k=1}^n [l_k(\theta) - q_k(\theta)] \\ & = q(\theta) + d(\theta) \end{align}$$

그리고 위 식의 $d(\theta) = d$ term을 estimate하는 것을 목표로 하며 최종적으로 log likelihood를 $\hat l_m = \sum_{k=1}^n q_k (\theta) + \hat d_m$와 같이 estimate 한다.

($l, q$의 차이로 $d$를 계산할 수 있다. 그리고 이것을 활용하여 $\eta_i$를 정의한다. $m$개 sampling한 eta의 sample mean을 통해 $\hat d_m$을 추정할 수 있으며, 최종적으로 이것을 활용하여 log likelihood를 estimate하는 것이다. 또한 $d$를 줄이면 $q$가 log likelihood에 approximate하게 되므로 make sense 하다.)

이는 결과적으로 서로 유사한 $(x_k, y_k)$ 데이터 셋은 유사한 $l_k(\theta)$값을 보이므로, 비슷한 데이터끼리 clustering하는 효과를 주며, $l_k(\theta)$의 second order Taylor expansion을 사용하여 구한 $q$는 각 cluster의 centroid 주변에 위치하게 된다. 그리고 위처럼 estimate한 log likelihood를 원래 scale로 변환하면 bias가 생기게 되어 이를 unbiased하게 만들기 위해 논문의 2.4 수식과 같이 약간의 term을 추가하여 다음과 같이 exponential scale로 변환한다 (call it as a bias-corrected estimator).

$$\hat p_m (y | \theta, u) = \exp\Big(\hat l_m(\theta) - \hat \sigma^2_\eta (\theta)/2m\Big)$$

2. Delayed acceptance MH

Delayed acceptance MH의 목표 역시 parameter $\theta$에 대한 posterior distribution $\pi(\theta) = p(y | \theta)p(\theta)/p(y)$로 부터 sampling을 하는 것이다. Payne and Mallick이 제안한 DAMH 알고리즘은 다음과 같다.

● Algorithm

1. First stage에서는 current parameter $\theta_c = \theta^{(j)}$와 첫번째 proposal distribution $q_1$로부터 $\theta^\prime \sim q_1(\theta | \theta_c)$을 propose하여 아래의 probability를 계산한다. 이 때 위에서 estimate한 likelihood function $\hat p$를 사용한다 (but without using control variates for $\hat l_m$).

$$\alpha_1(\theta_c \rightarrow \theta^\prime) = \text{min} \Big\{1, \frac{\hat p_m(y| \theta^\prime, u) p(\theta^\prime) / q_1(\theta^\prime | \theta_c)} {\hat p_m (y | \theta_c,u) p(\theta_c) / q_1(\theta_c | \theta^\prime)}\Big\}$$

2. 위의 First stage를 거친 후 $\theta_p$를 $\alpha_1(\theta_c, \theta^\prime)$의 확률로 $\theta^\prime$로 propose하고 또는 $1-\alpha_1(\theta_c, \theta^\prime)$의 확률로 $\theta_c$로 propose 한다. 그 이후 아래 확률을 계산하여 최종적으로 $\theta^{(j+1)}$를 update한다. (First stage에서 reject가 일어난 $\theta_p = \theta_c$인 경우 확률이 1이되므로 $\theta^{(j+1)} = \theta^{(j)}$가 된다.)

$$\alpha_2(\theta_c \rightarrow \theta_p) = \text{min} \Big\{1, \frac{p(y | \theta_p) p(\theta_p) / q_2(\theta_p | \theta_c)}{p(y| \theta_c) p(\theta_c) / q_2(\theta_c | \theta_p)} \Big\}$$

Stage 2의 proposal $q_2$는 spike and slab(first stage의 $\alpha_1$을 weight로 함) density로 설정 가능하며, Payne and Mallick이 사용한 위 알고리즘의 first stage에서 estimated likelihood $\hat p_m$은 current state와 independent하기 때문에 standard MH theory를 따른다. 이 방법은 independent함을 사용하여 theoretical 하게 true posterior로의 수렴의 증명이 더 쉽다는 장점이 있다. 만약 이를 $\hat p_m(y| \cdot,u) = \hat p^{(\theta_c)}_m(y | \cdot, u)$와 같이 state dependent하게 만든다면 standard MH theory를 적용할 수는 없지만 parameter의 posterior distribution $\pi(\theta)$에 converge한다는 것이 입증되어 좋은 측면이 있다. 따라서 저자들은 후에 delayed acceptance block PMMH 알고리즘에서 이를 활용한다.

또한 위 알고리즘은 이들이 사용한 control variates 방법을 사용하지 않았는데, 이로 인하여 sample의 variance가 커지게 되어 결과적으로 estimated likelihood $\hat p_m$에 악영향을 끼치게 된다고 주장한다. 저자는 만약 control variates를 사용하여 likelihood를 estimate한다면 bias correction 과정 없이도 likelihood를 더 잘 estimate한다고 주장한다.

● 알고리즘 성능 측정

알고리즘의 성능을 비교하기 위하여 첫번째는 statistical efficiency가 있다. 이는 estimation이 더 작은 variance를 가지면 더 효율적인 것으로 여겨지며 estimate가 true 분포에 얼마나 근접한지를 의미한다. 두번째는 computational efficiency가 있다. 이는 주어진 iteration 동안 얼마나 효과적인 sample을 만들어내는지를 의미한다.

기존의 MH에 비해 DAMH는 less statistically efficient하다. 그 이유는 추가적인 수락 확률을 계산해야하기 때문이다. 두 번 거르는 작업을 하기 때문에 DAMH는 MH에 비해 새로운 sample을 accept할 확률이 상대적으로 낮아지며 따라서 더 강한 autocorrelation을 가질 확률이 높아진다. 이들은 Theorem 1을 통해 두번째 stage에서의 accpetance probability가 $\sigma^2_R$에 따라 수치적으로 변하는 양상을 설명한다.

반면 statistical and computational efficiency의 균형적인 측면에서는 DAMH가 MH에 비해 뛰어나다. DAMH는 첫번째 stage를 통해 수락 확률이 낮은 sample은 거르기 때문에 computational efficiency가 조금 더 나은측면이 있으며 sample의 variance를 줄일 수록 statistical efficiency 역시 증가하기 때문이다. 구체적인 성능 측정 지표는 논문의 equation 3.7 ~ 3.9에 나와있다.

3. Real data application

저자들은 자신의 control variates를 사용한 방법과 Payne and Mallick이 제안한 방법을 real data application을 통해 비교하였다. 각 방법을 사용해 스웨덴 기업의 1991~2008년도 기간동안의 534,717개 데이터를 covariate으로 하여 파산 확률을 모델링하였다. 구체적으로 아래와 같은 logistic regression model을 사용하고 beta의 prior distribution으로는 multivariate Gaussian distribution을 사용하였다.

$$p(y_k | x_k,\beta) = \frac{\exp(x_k^\top \beta)}{1+\exp(x_k^\top\beta)}$$

각 알고리즘의 수행 결과로 Effective Draws (ED) $= \frac{N}{IF \times t}$를 계산하여 이것을 각각의 분자로 하고, standard MH의 ED를 분모로 하여 Relative ED (RED) 값을 비교해 성능을 측정하였다. 따라서 table의 RED가 높을수록 성능이 더 뛰어나다고 할 수 있다. 이들의 알고리즘을 사용한 결과인 Table 1 과 Payne and Mallick의 방법을 사용한 Table 2를 비교해보면 Table 1의 RED값이 더 높으며 특히 Payne and Mallick의 방법은 매우 적은 수의 subsample을 사용한 경우 성능이 좋지 않았다. Table 1의 결과를 보면 더 많은 수의 subsample을 사용하면 두번째 stage의 acceptance probability 또한 높아져 알고리즘이 더 효율적인 결과를 낸다는 것을 확인할 수있다.

Reference:

Quiroz, M., Tran, M. N., Villani, M., & Kohn, R. (2018). Speeding up MCMC by delayed acceptance and data subsampling. Journal of Computational and Graphical Statistics, 27(1), 12-22.

Link:https://www.tandfonline.com/doi/full/10.1080/10618600.2017.1307117

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.

1. Abstract

Spatio-temporal point process(STPP)란 시간, 공간에 따른 사건의 stochastic collection이다. 계산 복잡도를 고려하여, 기존의 시공간 통계학에서는 시간과 공간(time and space)의 conditional independence를 활용하여 distribution을 따로 구해 STPP의 solution을 도출하였다 (시간을 modeling하기 위해 Poisson process, Hawkes process 등이 사용되며 공간을 modeling하기 위해 kernel density estimator가 사용된다).

그러나 이들은 시공간을 통합하는 joint distribution을 구하지 못하는 것이 제한된 성능을 발휘할 수 밖에 없다고 생각한다 (예를 들어 지진의 발생은 이전 지진과 가까운 지역 및 시간대에 발생한다). 따라서 이를 극복하기 위해 STPP의 parameteriazaion framework를 diffusion model을 사용하여 제안한다.

Diffusion model의 input으로 spatial, temporal domain의 embedding 값을 구하여 사용하였으며, diffusion model을 통해 최종적으로 그들이 얻고자 하는 것은 hidden representation이 given 되었을 때 origial time and space를 estimate한 time and space이다.

2. Difficulties to derive spatio-temporal joint distributions

시공간의 joint distribution을 구하는 것은 아래와 같은 두가지의 이유로 쉽지만은 않다.

1. 시공간 joint distribution for STPPs는 sample space가 매우 거대하며, 이는 intractable하다. 따라서 우선 temporal density $p^*(t)$를 fitting한 후에 conditional density $p^*(s\mid t)$를 순차적으로 fitting하는 방법이 있으나, $p^*(s \mid t)$는 특정 model structure (KDEs, CNFs 등)에 제한된다.
2. 사건의 발생은 시공간 사이의 복잡한 correlation과 연관이 있다. 다양한 분야에서 발생한 사건은 그 분야에 따라 distinct한 spatio-temporal dependencies를 가진다. 따라서 대부분의 문제에 들어맞는 근본적인 시공간 dependence를 capture하는 방법이 필요하다.

따라서 이들은 새로운 parameterization framework인 Spatio-Temporal Diffusion Point Processes (DSTPP)를 제안한다. 간단히 아이디어를 설명하자면 diffusion 방법의 denoising probabilistic modeling을 활용하여 original complex distribution을 Markov chain of multiple steps로 decompose한다.

3. Spatio-temporal diffusion point processes

이들이 제안한 framework는 위와 같다. 크게 2개의 key module, spatio-temporal self-attention encoder와 spatio-temporal diffusion model로 구성되어 있다.

3.1 Spatio-temporal Encoder

이 단계는 spatio-temporal data를 embedding하고, hidden representation을 추출하는 과정이다. 시공간 데이터를 input하면 embedding layer와 attention layer를 거쳐 hidden representation을 반환한다. 저자들은 이 과정을 Encoding으로 생각하여  Encoder라고 표현한 듯 하다.

우선 시간을 embedding하는 방식은 아래와 같다. $M$은 embedding dimension이며 따라서 $e_t$가 temporal embedding을 나타낸다.

$$[e_t]_j = \begin{cases} \cos(t/10000^{\frac{j-1}{M}}) \text{ if j is odd} \\ \sin(t/10000^{\frac{j-1}{M}}) \text{ if j is even} \end{cases}$$

공간은 다음과 같이 embedding한다. $W_e$는 학습가능한 parameter를 포함하며 공간이 continuous domain인 경우 (like 좌표 데이터)와 discrete domain인 경우 (like lattice data)에 모두 사용 가능하다.

$$e_s = W_e s$$

또한 위 두개의 embedding 결과를 더하여 spatio-temporal embedding $e_{st}$를 정의할 수 있다.

결과적으로, 시공간 데이터 $S = \{(t_i, s_i)\}_{i=1}^L$에 대해 시공간 embedding $E_{st} = \{e_{st,1},\dots,e_{st,L}\}$, 시간 embedding $E_s = \{e_{t,1}, \dots, e_{t,L}\}$, 그리고 공간 embedding $E_s = \{e_{s,1}, \dots, e_{s,L}\}$을 얻을 수 있다.

그 다음 목표는 이 3가지의 embedding을 활용하여 서로다른 시공간마다의 characteristic을 capture하는 것이다. 따라서 논문의 (7), (8)번 수식을 통해 embedding vector들의 연산을 하고 최종적으로 spatial, temporal, and spatio-temporal 각각에 대한 hidden representation에 대해 얻을 수 있게 된다.

3.2 Spatio-temporal Diffusion and Denoising Processes

위에서 hidden representation, $h_{i-1}$을 얻은 그 다음의 목표는 hidden representation을 conditioning하여 미래 event에 대한 spatio-temporal joint distribution의 model을 학습하는 것이다. 그리고 여기에 diffusion model을 활용하였다.

이들은 space와 time에 각각 noise를 더하고 이를 복원하는 식으로 diffusion process와 reverse process를 진행하였다. 각 event $x_i =(\tau_i, s_i), \text{ (where } \tau_i \text{ denotes the time interval since the last event) }$에 대해 diffusion process modeling을 위해 K번의 diffusion process를 거친 Markov process $(x_i^0, x_i^1, \dots, x_i^K)$를 적용하였다. 각 diffusion step에 Guassian noise를 더해가는 과정을 다음과 같이 쓸 수 있다.

$$q_{st}(x_i^k \mid x_i^{k-1}) = (q(\tau_i^k \mid \tau_i^{k-1}),q(s_i^k \mid s_i^{k-1})), \\q(x^k \mid x^{k-1}) = N(x_k ; \sqrt{1-\beta_k}x^k, \beta_k \mathbf{I})$$

반면 noise를 제거해가며 원본 data $x_i = (\tau_i, s_i)$를 복원하는 과정, reverse process는 hidden representation을 condition으로 아래의 분포를 얻는 과정이다.

$$p_\theta(x_i^{k-1} \mid x_i^k, h_{i-1}) = p_\theta(\tau_i^{k-1} \mid \tau_i^k, s_i^k, h_{i-1})p_\theta(s_i^{k-1} \mid \tau_i^k, s_i^k, h_{i-1})$$

이전 event까지의 hidden representation을 given하여 다음 event의 시간과 공간의 joint distribution을 서로 conditionally independent하게 만들 수 있다는 의미인 듯 하다.

따라서 위 방법을 사용하여 diffusion model을 train하고, trained model을 사용하여 real data $(\tau_i^0, s_i^0)$에 대한 prediction값을 얻을 수 있다. 또한 이 prediction에서 중요한 것은 co-attention denoising network를 사용했다는 점이다. 이는 spatial과 temporal domain의 상호의존성을 capture하기 위함이며 3.4에 자세히 나와있다.

4. Experiments

이들은 real dataset과 synthetic dataset에서 다양한 실험을 진행하였다. 그 중 몇가지의 과정과 결과만 살펴보자. 이들은 Earthquakes, COVID-19 등의 real datasetrhk HawkewGMM의 synthetic dataset을 사용하였으며 이들이 구한 model vs 기존의 baseline model에서의 모델 성능을 두가지 측면으로 비교하였다. 첫번째는 space와 time각각의 negative log-likelihood를 구하여 비교하는 것이고, 두번째는 time은 RMSE, spatial location은 Euclidean distance를 계산하여 비교한 것이다.

결과를 살펴보면 이들의 model DSTPP는 모든 dataset에서 가장 뛰어나거나 두번째로 뛰어난 성능을 보였다. 기존의 NSTPP와 같은 conditional distribution을 사용하여 시공간의 distribution을 구하는 방식은 spatial 측면에서는 시간과 공간의 independence를 가정하고 modeling한 결과보다는 뛰어났지만, temporal 측면에서는 temporal modeling만 사용한 다른 방식들에 비해 아주 뛰어난 성능을 보이지는 않았다. 이러한 방법도 시공간의 joint distribution을 구한 것은 동일하지만 temporal domain은 잘 학습하지 못함을 보여주는 결과였다.

References:

Yuan, Y., Ding, J., Shao, C., Jin, D., & Li, Y. (2023, August). Spatio-temporal diffusion point processes. In Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining (pp. 3173-3184).

0. Overview

공간통계학에서 우리는 통계적 Modeling을 위해 보통 stationary 가정을 한다. 이는 간단히 말해서 특정 지점사이의 거리가 같다면, 공분산이 같음을 의미한다 (분산이 특정 지점 사이의 거리차이에 depend). 그러나 단순히 홍대입구역부터 신촌역까지의 거리와, 홍대입구역부터 합정역 까지의 거리의 같다는 사실이 각 지점의 특성이 같다는 것을 의미하지는 않는다. 예를 들어 강수량과 같은 수치는 단순한 거리 차이보다 각 지점의 온도나 고도 등에 더 많은 영향을 받는다.

따라서 이러한 공간적 특성을 반영하기 위해 modeling시에 nonstationary 가정을 해야하지만 이는 분석이 쉽지 않다. 이들은 Convolutional neural networks (CNN)을 사용하여 공분산을 modeling하였는데 이에 대해 살펴보도록 하겠다. 또한 CNN을 사용하기 위해 좌표 data를 어떻게 CNN의 input으로 변환하였는지를 중점적으로 살펴보도록 하겠다.

1. Methods

기존에 존재하는 nonstationary covariance modeling 방법은 1) spatial deformation method, 2) process convolution approach 등이다. 구체적인 설명은 생략하고 2)에서는 matern covariance의 closed form을 구할 수 있는데, 지역의 수가 많아질수록 계산해야하는 matrix의 크기가 커지며 따라서 추정하기 힘들다는 단점이 있다. 따라서 구역(subregion)을 나누어 local stationary를 가정하는 방법들이 후에 제안되었는데 이러한 방법은 subregion을 나누는 것이 주관적이라는 단점이 존재한다.

이들은 spatial data를 incomplete image로 여기고(image는 특정 image size에서 모든 점들이 관측되어 RGB의 값을 가지지만, spatial data는 관측이 되지 않는 지점이 존재하므로 이를 관측되지 않은 RGB로 생각) CNN을 사용하여 주관적으로 subregion을 나누는 기존의 방법을 data-driven 방법을 제안함으로서 (조금 더 객관적) 개선하고자 한다. 

2. CNN model

2.1  CNN architecture

CNN model을 사용한 이유는 위에서 언급한 것과 같이 stationary 와 nonstationary region을 구분하기 위함이다. 따라서 100 by 100 크기의 input에 따라 구역이 nonstationary할 확률을 output한다. 따라서 분류문제로 생각하여 확률이 0.5보다 크거나 같으면 저자는 nonstationary 구역으로 생각하고 그렇지 않으면 stationary 구역으로 생각하여 모델을 학습하고 사용한다.

구체적으로 이들이 사용한 CNN 구조는 위와 같다. 우선 100 by 100 size의 input data를 3 by 3의 kernel 32개를 사용하여 feature를 extract하는 convolution layer를 거친다. 그 후 flatten layer로 이를 vector로 펼치고 마지막으로 FC layer에서 두 단계의 non-linear transformation을 거쳐 최종적으로 nonstationary region일 확률을 나타낸다.

2.2 Data pre-processing

그러나 문제점이 존재하는데 CNN을 사용하기 위해서 100 by 100 size의 input 형태로 data를 변환해야하지만 현실 세계의 많은 데이터는 100 by 100 grid에 맞게 균일하게 위치하지 않는다는 것이다. 따라서 이들은 3단계의 전처리 과정을 거쳐 100 by 100 size의 input형태를 만들어내며 구체적인 설명은 다음과 같다. 

• Splitting: divide observed region $[0,1]^2$ to 100 by 100 subregions (denote each subregion as $D_{i,j}$)
• Averaging: $D_{i,j}$에 속하는 observation의 개수를 $N_{i,j}$, $D_{i,j}$에 속하는 observation들의 평균을 $\bar{Z}(i,j)$라고 하여 아래와 같이 계산한다.

$$\bar{Z}(i,j) = \begin{cases} \frac{1}{N_{i,j}} \sum_{k:s_k \in D_{i,j}} Z(s_k), \quad N_{i,j}>0 \\ \bar{Z}(u,v), \ \text{where }(u,v) = argmin_{u,v}\{(i-u)^2 + (j-v)^2 : N_{u,v}>0\}, \quad N_{i,j}=0 \end{cases}$$

• Scaling: 위에서 얻은 $\bar{Z}(i,j) \in [0,1]$의 minimum을 m, maximum을 M이라고 하자. 아래와 같이 scaling을 한다.

$$\tilde{Z}(i,j) = \begin{cases} (\bar{Z}(i,j) - m) / (M-m), \ \ m \neq M; \\ 0.5, \quad \quad m=M \end{cases}$$

위 과정을 요약하자면, $[0,1]^2$의 좌표를 100 by 100의 subregion으로 나눈 후, 각 subregion에 속하는 점들의 GRF값의 평균을 취한 다음 minmax scaling을 한다. 결과적으로 우리가 원하는 100 by 100 size의 input 형태를 갖추게 된다. (그런데 location 좌표 $s_k$를 subregion $D_{i,j}$에 어떻게 mapping하는 지는 나와있지 않다)

2.3 Subregion selection

이제 남은 것은 train한 CNN과 전처리한 data를 사용하여 stationary한 지역과 nonstationary한 지역을 나누는 것이다. 이를 위해서 subregion을 나누어야 한다 (input size를 100 by 100으로 나누기 위해 2.2에서 정의한 subregion과 다름). 이들이 제안한 방법을 정리하면 다음과 같다.

1.  좌표 $S = \{s_1,\dots,s_n\}$에서 K개의 random한 point $A = \{a_1, \dots, a_K\}$를 선택한다.
2. 그 후 distance vector $E = \{ \lVert s_i - a_k \rVert^2: a_k \in A\}$를 계산한다. distance vector는 특정 point $s_i$에서 random하게 선택된 K개의 point $a_k$까지의 거리이며 따라서 성분이 K개이다.
3. 특정 점 $s_i$가 $a_j$와 거리가 가장 최소라면, subregion $R_j$에 $(s_i, Z(s_i))$를 할당한다.
4. 각 subregion$\{R_k\}_1^{K}$마다 2.1의 CNN 구조를 사용하여 nonstationary probability를 구한다 (let this $p_k, k=1,\cdots,K$).
5. $P = \sum_{k=1}^K p_k$를 구한 후 iteration을 반복하여 최소가 되는 $P$를 찾고 그때의 subregion에 속하는 점들의 centre를 반환한다 (각 iteration마다 선택되는 점의 집합 $A$가 다르므로 $P$도 달라진다).

요약하자면 몇 개의 subregion으로 나눌 것인지만 선택하면, 우리가 가진 데이터들이 어느 subregion에 속하는지 알 수 있는 알고리즘이다. 다만 최적의 subregion수 K를 선택하는 방법이 없다는 것은 조금 아쉽다.

3. Simulation

이들은 두 종류의 simulation을 하였다. 간략히 설명하자면 첫번째 실험은 stationary or nonstationary여부를 CNN이 잘 구분하는가에 관한 실험이다. 과정은 우선 stationary, nonstationary dataset을 각각 16000개 generate 한 후 80%의 sample을 활용해 CNN을 train한다. 구체적인 data generation과정은 생략하고, 남은 20%의 데이터를 활용해 test한 결과는 아래와 같다.

간혹 stationary data인 것을 nonstationary로 판단하는 경우(and 그 반대)가 있었지만 대체로 정확한 classification을 한 것을 알 수 있다.

두번째 실험은 $[0,1]^2$ domain에서 4개의 기준 점(anchor locations)을 기준으로 하여 4가지의 서로 다른 nonstationary parameter를 가지는 subregion에서 10,000개의 nonstationary data를 생성한다. 그 후 2&3개의 subregion을 가지도록 CNN을 학습한다. 여기서 실제로 나눈 subregion과 비슷하게 구역이 나뉘는지 check할 수 있으며, ExaGeoStat package를 사용하여 nonstationary data의 parameter또한 추정할 수 있는데, CNN model을 사용하여 추정한 parameter의 MSE vs. 기존 방식으로 subregion을 나눈 후 추정한 parameter의 MSE를 비교하였다.

추정 결과를 간단히 살펴보자면 CNN으로 subregion을 나눈 후 parameter를 추정하면 user-defined subregion보다 실제 parameter를 더 잘 추정하는 것을 알 수 있다.

4. 요약

이 논문은 공간통계학에서 stationary 가정이 위배되는 많은 상황이 있다는 것에서 출발하여, 기존의 local stationary 방법이 주관적인 선택에 의해 subregion을 만드는 것에 문제를 제기하고, CNN을 통해 subregion을 만드는 방법을 개선하고자 하였다. Stationary 가정을 local하게 하되, 객관적인 방법으로 subregion을 만든 후 subregion내의 point 사이의 local stationary를 가정하여 분석하는 방안을 제시한 것이다.

References:

Nag, P., Hong, Y., Abdulah, S., Qadir, G. A., Genton, M. G., & Sun, Y. (2023). Efficient Large-scale Nonstationary Spatial Covariance Function Estimation Using Convolutional Neural Networks. arXiv preprint arXiv:2306.11487.

0. Preliminary

ETAS point process

시작에 앞서 논문을 이해하기 위해 ETAS point process에 대한 사전지식이 필요하다. Epidcmie-Type Aftershock Sequence (ETAS) model이란 날짜, 시간, 위도, 경도, 지진의 진도 등을 사용하여 지진의 발생을 modeling하는 과정이다.

수식적인 설명을 덧붙이면 이는 conditional intensity $\lambda(t, x, y, m\mid H_t)$를 사용한 point process 이다.

이 intensity function의 notation들을 조금 더 자세히 살펴보자면, $(x, y)$는 2차원 좌표이며 $m, t$는 진도와 시간이다. 또한 $H_t$는 $t$ 시점까지의 process들의 history이며 $\Big(H_t = (t_i, x_i, y_i, m_i) : t_i \le t \Big)$, 따라서 intensity function은 t시점 이전까지의 event들이 given되었을 때, 시간 $t$, 좌표 $(x,y)$, 진도 $m$인 이벤트가 발생할 expected rate라고 해석할 수 있다.

또한 이는 Hawkes' point process의 special case인데, conditional intensity function을 아래와 같이 decompose 할 수 있다.

$$\lambda_{\beta, \Theta}(t, x,y,m\mid H_t) = \nu_\beta(m)\lambda_\Theta(t,x,y \mid H_t)$$

여기서 $\nu_\beta (m)$은 magnitude m 의 pdf with unknown parameter $\beta >0$이다. 오른쪽 $\lambda$ 함수는 다음과 같이 decompose가 가능하다.

$$\lambda_\Theta(t,x,y\mid H_t) = \tilde{u}(x,y) + \sum_{t_i < t} \kappa_{A,\alpha}(m_i) g_{c,p}(t-t_i) f_{D,\gamma,q}(x-x_i, y-y_i;m_i)$$

$\tilde{u} (x,y)$: $\mu >0$을 모수로 가지는 background seismicity (지진활동도) rate를 나타낸다.

$\kappa_{A,\alpha}(m_i)$: 진도가 $m_i$로 발생하는 지진의 기대 발생 횟수이다.

$g_{c,p}(t-t_i)$: 지진 발생 시간에 대한 probability density function이다.

$f_{D, \gamma, q}(x-x_i, y-y_i ; m_i)$: 지진이 발생하는 좌표에 대한 probability density function이다. 이때 진도는 $m_i$이다.

Interpretable Deep Generative Spatio-Temporal Point Processes

1. Goal

위 논문에서는 Neural Embedding Spatio-Temporal (NEST) point process를 제안한다. 위에서 언급한 지진의 모델링에 많이 사용하던 ETAS는 kernel function이 location에 상관 없이 homogeneous한 성질이 있어서 model fitting이 용이하였지만, 이러한 고전적인 모델들은 conditional intensity의 parameter 값에 의존할 수 밖에 없다는 단점 또한 존재한다.

따라서 NEST model은 복잡한 spatial dependence를 더 잘 capture하고, 해석이 가능하며, computation에도 장점이 있다.

(intensity function이란 시간, 공간에 따른 사건 발생률의 함수이며, conditional intensity function은 과거 사건의 발생을 고려한다.)

2. Method

이들은 continuous-time and continuous space point process model인 NEST model을 제안한다. NEST model에서는 intensity function으로 neural network를 사용하게 되는데, parameter에 의존하는 기존의 방법들과는 달리 이들은 neural network를 사용함으로 model의 interpretability를 유지하며 flexibility를 줄 수 있다고 주장한다.

구체적으로 이들은 다음과 같은 conditional intensity function을 제안한다.

$$\lambda^*(t,s) = \lambda_0 + \sum_{j:t_j<t}\nu (t,t_j,s,s_j)$$

$\nu$는 과거 사건 $\mathcal{H}_t$의 영향을 capture하는 kernel function이며, 여기서는 standard Gaussian diffusion kernel 형태의 커널을 사용했다. (사건 간의 시공간적 dependency를 결정하는 함수)

kernel function의 구체적인 형태는 다음과 같다.

$$\nu(t,t^\prime, s, s^\prime) = \sum_{k=1}^K \phi_{s^\prime}^{(k)} \cdot g(t,t^\prime, s, s^\prime \mid \Sigma_{s^\prime}^{(k)} , \mu_{s^\prime} ^{(k)}), \quad \forall t^\prime < t, s \in \mathcal{S}$$

term들을 조금 더 자세히 살펴보자면,

$\{\mu_{s^\prime}^{(k)}, \Sigma_{s^\prime}^{(k)}\}$:  mean and covariance matrix parameters

$K$:  Gaussian mixture의 component 수를 결정하는 hyper-parameter

$\phi_{s^\prime}^{(k)}:\mathcal{S} \rightarrow \mathbb{R}$:  $\sum_{k=1}^K \phi_{s^\prime}^{(k)}=1, \forall s^\prime \in \mathcal{S}$를 만족하는 k-th Gaussian component의 weight

$g(\cdot)$ 함수의 구체적인 형태는 다음과 같다.

$$g(t, t^\prime, s, s^\prime) = \frac{Ce^{-\beta(t-t^\prime)}}{2\pi \sqrt{\lvert \Sigma_{s^\prime} \rvert}(t-t^\prime)}\cdot \exp \Big\{ - \frac{(s-s^\prime - \mu_{s^\prime})^\top \Sigma_{s^\prime}^{-1}(s - s^\prime - \mu_{s^\prime})}{2(t-t^\prime)} \Big\}$$

$\beta>0$가 커질수록 위 함수의 값이 빠르게 감소하므로, 이는 temporal decay rate를 조절하는 parameter이며, $C>0$은 constant로 값의 magnitude를 조절한다.

$\mu_s = [\mu_x(s), \mu_y(s)]^\top, \Sigma_s$는 diffusion kernel의 mean and covariance parameters를 의미하며, Ganssian component들의 shift, rotation, shape을 결정하는 parameter이다.($\Sigma$는 positive semi-definite matrix이다).

3. Conclusion

결과적으로 저자들은 자신들의 방법 (NEST)가 기존의 ETAS 방법을 사용했을 때보다 spatial dependency를 더 잘 capture한다고 주장한다. 또한 model parameter 추정에서 neural network를 사용한 것 같은데, model의 기조는 Gaussian diffusion kernel이므로 parameter값에 대한 해석이 가능하기 때문에 interpretable이라는 단어를 논문 제목에 붙인 것 같다. 구체적인 실험 결과 및 수식은 본문의 Chapter 4와 Appendix에서 확인 가능하다.

References:

Zhu, S., Li, S., Peng, Z., & Xie, Y. (2020). Interpretable deep generative spatio-temporal point processes. In AI for Earth Sciences Workshop at the 34th Conference on Neural Information Processing Systems (NeurIPS 2020), online, Dec(pp. 6-12).

Ogata, Y. (1998). Space-time point-process models for earthquake occurrences. Annals of the Institute of Statistical Mathematics, 50, 379-402.