060. Parameter Estimation of ERGM

Parameter Estimation of ERGM

Current Methods for ERGM

1. Approximation-based Algorithm:
   MCMC 샘플들로 likelihood function 을 최대화. typically used in finding MLE.

   • Maximum Pseudo-likelihood Estimation (MPLE, 1974).
   • Markov Chain Monte Carlo Maximum Likelihood Estimation (MCMCMLE, 1994).
   • Markov Chain Monte Carlo Stochastic Approximation (MCMCSA, 2002).
     • 이 위까지로는 model degeneracy issue 를 해결하지 못했음. 아래의 VTSAMCMC 가서야 해결함.
   • Varying Truncation Stochastic Approximation MCMC Method with Trajectory Averaging (VTSAMCMC, 2013).

2. Auxiliary Variable Markov Chain Monte Carlo (MCMC) Algorithm:
   normalizing constant 를 근사하거나, normalizing constant ratio $\frac{\kappa (\theta )}{\kappa ({\theta }^{\prime })}$ 를 cancel 하기 위해 auxiliary variable 들을 도입 (introduce). 베이지안 추론을 위해 사용. Bayesian inference.

   • The Exchange Algorithm.
   • Auxiliary Variable Metropolis-Hasting Algorithm (AVMH).
   • Adaptive Exchange Monte Carlo Algorithm (AEXMC).

Approximation-based Algorithm

Maximum Pseudo Likelihood Estimation (MPLE)

$A$ 내부의 성분 간의 의존성 (high-order dependence wihin network) 을 무시한 채로, conditional independence 를 가정하여 조건부 likelihood 함수들의 series를 곱하는 것으로 pseudo-likelihood 함수를 근사함. ERGMS 의 조건부이며 pseudo-likelihood 는 아래와 같음. 이렇게 근사한 pseudo-likelihood 로 maximum likelihood estimator 를 approximate. 물론 이건 어디까지나 가장 간단한 방법으로서 제시되었을 뿐이고, 이 방법론은 의존성 구조를 완전히 무시했기 때문에 퍼포먼스가 구림.

$$\begin{array}{rlrl}& logit\left\{P_{\theta }\right(X_{ij}=1|X_{ij}^c=x_{ij}^c)\}& & ={\theta }^{\prime }g(x_{ij}^c)\\ ⟺& \log PL(\theta ,x)& & =\sum _{ij}{\theta }^{\prime }g(x_{ij}^c)\cdot x_{ij}-\sum _{ij}\log \left\{1+{\theta }^{\prime }g(x_{ij}^c)\right\}\end{array}$$ $$\begin{array}{rlrl}& logit\left\{P_{\theta }\right(X_{ij}=1\left|y_{ij}^c\right)\}& & =\underset{\text{log of unnormalized density with yijc}}{\underbrace{\sum _{i=1}^P{\theta }_i^{\prime }{S}_i(y_{ij}^c)}}\\ ⟺& \log PL(\theta ,y)& & =\underset{\text{not a proper likelihood}}{\underbrace{\sum _{ij}{\theta }^{\prime }S(y_{ij}^c)\cdot y_{ij}-\sum _{ij}\log \left\{1+{\theta }^{\prime }S(y_{ij}^c)\right\}}}\end{array}$$

• advantage: simple one / don’t need to generate auxiliary networks
• disadvantage: parameter estimates is not good because MPLE ignores high-order depedence

MCMC MLE

1. with initial guess of parameter → generate auxiliary networks
2. with auxiliary networks → update parameter estimates

패러미터 ${\theta }^{(0)}$ 에 대해 이하와 같은 initial guess 가진다고 하자.

$$\begin{gathered} l(\theta )={\theta }^{\prime }S(y)-\log k(\theta ) \\ l({\theta }^{(0)})={{\theta }^{(0)}}^{\prime }S(y)-\log k({\theta }^{(0)}) \end{gathered}$$ $$l(\theta )-l({\theta }^{(0)})=(\theta -{\theta }^{(0)}{)}^{t}S(y)-\log \frac{\kappa (\theta )}{\kappa ({\theta }^{(0)}}$$

distribution $f(x|{\theta }^{(0)})$ 에서 생산된 MC 샘플들로부터 $\kappa (\theta )$ 근사. 이때 ${\theta }^{(0)}$ 는 $\theta$ 의 initial 측정값. MCMC 시뮬레이션을 통해 $f(x_{obs}|{\theta }^{(0)}$ 으로부터 랜덤샘플 $x_1,\cdots ,x_m$ 을 draw. 이 경우 log-ratio likelihood 는 $m\to \infty$ 에서 이하와 같이 근사 가능.

$$I(\theta )-I({\theta }^{(0)})=(\theta -{\theta }^{(0)}{)}^{t}g(x_{\mathrm{obs}})-\log \left[\frac{1}{m}\sum _{i=1}^m\exp \left\{(\theta -{\theta }^{(0)}{)}^{t}g(x_i)\right\}\right]$$

따라서 $\theta$ 를 경유해서 $I(\theta )-I({\theta }^{(0)})$ 를 극대화하는 작업을 통해 결국 ${\hat{\theta }}_{MLE}$ 를 근사할 수 있게 된다.

MCMLE 의 퍼포먼스는 ${\theta }^{(0)}$ 의 선택에 좌우된다. ${\theta }^{(0)}$ 이 MLE 의 참값의 attraction region 에 존재하지 않는다면, 해당 방법론은 suboptimal 해로 수렴해버리거나, 최악으로는 converge 실패 가능성도 있음.

MCMC Stochastic Approximation

exponential family 에 대한 이론을 통해, 우리는 ERGMs을 maximize 하는 것은 $E_{\hat{\theta }}\{g(X)\}=g(x_{obs})$ 라는 system 을 푸는 것과 동일하다는 것을 알 수 있다. 이 system 의 성립 여부는 이하와 iff ($⟺$).

1. $\hat{\theta }={\hat{\theta }}_{MLE}$
2. Independent network generation
3. Parameter estimation update with stochastic approximation

이 방법론은 independent 네트워크 샘플을 생산하는데에는 비효율적이다. 각 샘플 $x_{k+1}$을 생산하기 위해 요구되는 이터레이션 스텝의 갯수는 order of $100n^2$. 이때 $n$은 node 의 갯수.

Auxiliary Variable MCMC-based Approaches

Exchange Algorithm

normalizing constant ratio $\frac{\kappa (\theta )}{\kappa ({\theta }^{\prime })}$ 따위의 auxiliary 변수로 분포 $f(x|\theta )$ 를 augment. 이는 시뮬레이션 중에 canceled.

1. 후보 point ${\theta }^{\prime }\sim q({\theta }^{\prime }|\theta ,x)$ 를 생산
2. perfect sampler 사용해서 auxiliary 변수 $y\sim f(y|{\theta }^{\prime })$ 를 생산
3. ${\theta }^{\prime }$를 with probability $1\wedge r(\theta ,{\theta }^{\prime }|x)$로 채택. 이때

$$r(\theta ,{\theta }^{\prime }|x)=\frac{\pi ({\theta }^{\prime })}{\pi (\theta )}\cdot \frac{f(x|{\theta }^{\prime })}{f(x|\theta )}\cdot \frac{f(y|\theta )}{f(y|{\theta }^{\prime })}\cdot \frac{q(\theta |{\theta }^{\prime },x)}{q({\theta }^{\prime }|\theta ,x)}$$

exchange 알고리즘은 Ising 이나 autologistic 과 같은 일부 discrete 모델에는 잘 작동. 하지만 perfect sampling 이 적용되지 않는 많은 다른 모델에는 적용할 수 없다. 우리는 MCMC 샘플을 통해 auxiliary 변수 $y$를 생산할 수 있지만, 수렴 문제 부분에서 이론적인 흠결이 있음. 만약 MCMC 샘플의 mixing이 매우 느리다면 딴놈이 아니라 ${\theta }_0$ 가 매우 높은 확률로 채택되어버림. 그리고 이 MCMC 샘플의 mixing 이 느린 것 자체가 ERGMs 의 일반적인 특징임. 이래서 문제.

Monte Carlo MH Algorithm

MH 알고리즘의 MC 버전. 각 이터레이션에서 MCMH 알고리즘은 unknown normalizing constant ratio $\frac{\kappa (\theta )}{\kappa ({\theta }^{\prime })}$ 를 MC estimate와 importance 샘플링 접근법으로 대체한다.

exchange 알고리즘과 다르게, MCMH 알고리즘은 perfect sampler 요구조건을 빗겨간다. 따라서 perfect sampler 를 못 쓰는 다수의 통계문제에 대해서도 얘를 쓸 수 있음. 하지만 얘도 여전히 수렴 문제가 존재함. 얘의 경우에는 importance sampling estimator 가 유한한 숫자의 샘플로는 ratio의 참값 $\frac{\kappa (\theta )}{\kappa ({\theta }^{\prime })}$ 에 수렴하지 못할 수도 있음.

knitr::include_graphics('pics/network/06/01.png')

Varying Trunction Stochastic Approximation MCMC

ERGM 의 likelihood 함수는

$$\begin{array}{rlrlrl}& f(y|\theta )=\frac{1}{\kappa (\theta )}\exp \{{\theta }^{\prime }S(y)\},& & \theta =({\theta }_1,\cdots ,{\theta }_d{)}^{\prime },& & S(y)=(S_1(y),\cdots ,S_d(y){)}^{\prime }\end{array}$$

이때 $\kappa (\theta )$를 특정할 수 없다 (intractability)는 점과 모델 degeneracy 때문에, $\theta$를 정확하게 추정하는 것은 어렵다. 이 문제는 varying truncation stochastic approximation MCMC 를 사용하는 것으로 해결 가능.

Normalizing constant 는 $\kappa (\theta )={\sum }_{\text{all possible }y}\exp \{{\theta }^{\prime }S(y)\}$.

이로 인해 촉발되는 문제는?

1. MPLE: dyadic 독립이 가정되어 있다. 이건 observed 네트워크가 무엇이냐에 대해 과하게 의존.
2. MCMLE: ${\theta }^{(0)}$ 을 어떻게 고르느냐에 대해 과하게 의존. local 최적해로 수렴해버리거나 모델 degeneracy 로 수렴 자체가 실패하기도.

MCMC Stochastic Approximation

$h(\theta )=0$ 형의 system 을 풀자. 고전적인 SA 알고리즘은 이하와 같은 형태를 띈다. ($a_k$ , $h$, ${\omega }_k$ 에 대한) 적절한 조건 하에서, 이 알고리즘은 해로 수렴한다는 것을 실제로 보이는 것이 가능하다.

$$\begin{array}{rlrl}{\theta }_{k+1}& ={\theta }_k+a_k{Y}_k& & \\ & ={\theta }_k+a_k\{h({\theta }_k)+{\omega }_k\},& & k\ge 0\end{array}$$

• $Y=h(\theta )+\omega$ 는 noisy estimate
• $\omega$ 는 mean-zero noise.

${\theta }_{MLE}$ 를 찾는 것은 exponential family 안에서 $E_{\theta }\{S(Y)\}=S({\mathbf{y}}_{obs})$ 의 해를 찾는 것과 equivalent. 이는 이하의 단계를 거친다. 다만 이는 independent 네트워크 샘플을 생산하는데 있어 비효율적일 수밖에 없음. order of $100n^2$ 이라 연산을 엄청 먹으니까.

1. ${\mathbf{y}}^{k+1}\sim f(\mathbf{y}|{\theta }^{(k)}$를 샘플링 (independence 네트워크 생산)

각 arc 변수 $Y_{ij}$가 독립적으로 정해지는 랜덤 그래프에서 시작. 이 변수에는 0 혹은 1이 0.5 확률로 할당. 이 랜덤 그래프를 Gibbs 샘플러든 MH 알고리즘이든 써서 업데이트.

2. SA 를 통해 패러미터 estimate ${\theta }^{(k+1)}={\theta }^{(k)}-a_k{D}^{-1}\{U({\mathbf{y}}_{k+1},{\bar{\mathbf{y}}}_{k+1})-\mathbf{S}({\mathbf{y}}_{obs})\}$ 절차 따라서 업데이트. where

$$\begin{array}{rl}U({\mathbf{y}}_{k+1},{\bar{\mathbf{y}}}_{k+1})-\mathbf{S}({\mathbf{y}}_{obs})\}& =P({\bar{\mathbf{y}}}_{k+1}|{\mathbf{y}}_{k+1})\cdot \mathbf{S}({\bar{\mathbf{y}}}_{k+1})\{1-P({\bar{\mathbf{y}}}_{k+1}|{\mathbf{y}}_{k+1})\}\cdot \mathbf{S}({\bar{\mathbf{y}}}_{k+1})\\ {\bar{\mathbf{y}}}_{k+1}& =1-{\mathbf{y}}_{k+1}\end{array}$$

Model Degeneracy

$\theta$을 어떻게 정하느냐에 따라 모델은 그것의 확률을 (Complete (fully connected) 거나 empty (entirely unconnected) 네트워크와 같은) 1개 혹은 소수의 그래프에만 한정시켜서 부어버릴 수도 있다. (탐색 효율 쓰레기됨) 이 문제를 해결하기 위해선 탐색 전에 degeneracy 리전을 거의 포함하지 않는 패러미터 스페이스를 가지는 모델로 특정할 필요가 있다. 근데 이게 진짜 드럽게 어려움.

Varying Truncation SAMCMC for ERGMs

알고리즘 proceeds 전 Setup: for some constants $k_0>1$, $\eta \in (\frac{1}{2},1)$, $\xi \in (\frac{1}{2},\eta )$, $C_a>0$, $C_b>0$, set as below.

• tag: C1

$$\begin{array}{rlrl}& a_k=C_a{\left(\frac{k_0}{(k_0\vee k)}\right)}^{\eta },& & b_k=C_b{\left(\frac{k_0}{(k_0\vee k)}\right)}^{\xi }\end{array}$$

Next,

• tag: C2

$$\bigcup _{s\ge 0}{\mathcal{K}}_s=\Theta ,\text{where }{\mathcal{K}}_s\subset \text{int}({\mathcal{K}}_{s+1})$$

And also,

• $\mathcal{X}$: a space of social network
• $\mathcal{T}$: $\mathcal{X}\times \Theta \to {\mathcal{X}}_0\times {\mathcal{K}}_0$ (reinitialization mechanism)
• ${\sigma }_k$: 이터레이션 $k$ 까지 수행되는 reinitialization 의 숫자. (${\sigma }_0=0$)

---

Proceeds:

1. Gibbs 샘플러 사용해 $m$ sweeps 번 interate 해서 auxiliary network $y_{k+1}\sim f(y|{\theta }^{(k)})$ 생산

2. Set ${\theta }^{(k+\frac{1}{2})}={\theta }^{(k)}+a_k\{S(y_{k+1})-S(y_{obs})\}$

3. Set

$$\{\begin{array}{ll}{\sigma }_{k+1}={\sigma }_k\text{ and }\left(y_{k+1},{\theta }^{(k+1)}\right)=\left(y_{k+1},{\theta }^{(k+\frac{1}{2})}\right)& \Vert {\theta }^{(k+\frac{1}{2})}-{\theta }^{(k)}\Vert \le b_k,{\theta }^{(k+\frac{1}{2})}\in {\mathcal{K}}_{{\sigma }_k}\\ {\sigma }_{k+1}={\sigma }_k+1\text{ and }\left(y_{k+1},{\theta }^{(k+1)}\right)=\mathcal{T}\left(y_k,{\theta }^{(k)}\right)& o.w.\end{array}$$

Trajectory averaging estimator

trajectory averaging estimator ${\bar{\theta }}_n=\frac{1}{n}{\sum }_{k=1}^n{\theta }^{(k)}$ 에 의해 $\theta$ 를 estimate 가능.

실전에서는 estimate 의 variation 을 줄이기 위해 대신 $\theta$ estimate 에 $\bar{\theta }(n_0,n)=\frac{1}{n-n_0}\sum _{k=n_0+1}^n{\theta }^{(k)}$ 을 자주 사용함. 이 때 $n_0$ 는 burn-in 이터레이션의 숫자.

• Free parameters
  • Free parameters: $\{a_k\}$, $\{b_k\}$, $\{{\mathcal{K}}_s,s\le 0\}$, $m$
  • $k_0=100$, $\eta =0.65$, $\xi =\frac{0.5+\eta }{2}$.
  • $C_a$, $C_b$: adjusted for different examples
  • choose ${\mathcal{K}}_0$ to be around MPLE.
  • 이 article 에선 ${\mathcal{K}}_{s,1}=[-4(s+1),4(s+1)]$, ${\mathcal{K}}_{s,2}=\cdots ={\mathcal{K}}_{s,d}=[-2(s+1),2(s+1)]$ 로 설정.
  • $m=1$

Numerical Examples

Methods: MCMLE, SAA (Stochastic Approximation Algorithm), Varying truncation SAMCMC SAMCMC: independent 5 runs(each of 200,000 iterations) (m=1, Ca = 0.01, Cb = 1000, η, ξ, κs are defalut)

?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\?????????????????\

Conclusion

SAMCMC 알고리즘의 varying truncation 은 ERGM 에 적용 가능. 여기서는 estimator 가 degeneracy region 으로 향하는 걸 reinitialization 에 의해 방지. estimator 는 ERGM 의 MLE로 converge 하며, 이는 asymptotic 하게 $N$ 을 따르며 asymptotic 하게 efficient. 이는 normalizing constant 를 감당하고서라도 model degeneracy problem 를 해결할 수 있다는 점에서 유용.