Sansan Builders Blog

Sansanのものづくりを支えるメンバーの技術やデザイン、プロダクトマネジメントの情報を発信

【Dive into the Economics of Networks】vol.8 Double Metropolis-Hastings Sampler によって ERGM のパラメーターの事後分布を推定する

f:id:shiori_baba:20220316114139j:plain

Sansan R&D 研究員の小松です。前回までは、ネットワークをシミュレーションする際に発生する問題、およびその解決方法の1つについて議論してきました。

buildersbox.corp-sansan.com

なぜこれらが重要かというと、効用関数のパラメーターを推定するためにネットワークをシミュレーションする必要があるからです。

A naive Metropolis-Hastings sampler doesn't work

今考えている尤度関数は以下のようなものでした。


\begin{align}
\pi(g, \theta) =  \frac{\exp(Q(g, \theta))}{c(\mathcal{G} ,\theta)}
\end{align}

ただし、 Q はポテンシャル関数、 c(\mathcal{G} ,\theta) :=  \sum_{\omega \in \mathcal{G}} \exp(Q(\omega, \theta)) です。 この normalizing constant  c(\mathcal{G} ,\theta) があるために、尤度関数を直接計算することが困難なのでした。

Bayesian estimation の際も、事後分布を求めるときに normalizing constant が計算できないというケースがあります。  \theta の事前分布を  p(\theta) とおくと、事後分布  p(\theta |g)


\begin{align}
p(\theta |g) =  \frac{\pi(g, \theta) p(\theta)}{\int_{\Theta}\pi(g, \theta) p(\theta)d\theta}
\end{align}

と表せます。この事後分布の分母は一般に計算が困難です。そこで Metropolis-Hastings 法によって、この事後分布から  \theta をサンプリングすることを考えます。

Metropolis-Hastings 法の手順は以下の通りです。

  1. パラメーター  \theta の初期値  \theta_{0} を定める。
  2. 提案分布  q(\cdot| \theta) から  \theta^{\prime} をドローする。
  3. Acceptance ratio  \alpha(\theta, \theta^{\prime}) := \min\left( 1, \frac{ p(\theta^{\prime} |g) }{p(\theta |g)} \frac{q(\theta | \theta^{\prime})}{ q(\theta^{\prime} | \theta)} \right) を計算する。
  4. 一様分布  U[0, 1] から  u をドローする。
  5.  u \leq  \alpha(\theta, \theta^{\prime}) ならば、 \theta \leftarrow \theta^{\prime} とする。
  6. 2 に戻る。

Metropolis-Hastings 法の詳細はここでは触れません。例えば鎌谷 (2020) をご覧になると良いでしょう。

Metropolis-Hastings 法を用いるモチベーションは、normalizing constant  c(\mathcal{G} ,\theta) の計算を回避することです。 ただし上のアルゴリズムで計算される acceptance ratio は、


\begin{align}
\alpha(\theta, \theta^{\prime}) = \min\left( 1, \frac{ p(\theta^{\prime} |g) }{p(\theta |g)} \frac{q(\theta | \theta^{\prime})}{ q(\theta^{\prime} | \theta)} \right)
= \min\left( 1, \frac{ \exp(Q(g, \theta^{\prime})) }{ \exp(Q(g, \theta))  } \frac{ c(\mathcal{G} ,\theta) }{ c(\mathcal{G} ,\theta^{\prime}) } \frac{ p(\theta^{\prime}) }{ p(\theta) } \frac{ q(\theta | \theta^{\prime}) }{ q(\theta^{\prime} | \theta)} \right)
\end{align}

となり、normalizing constant が消えません。これは、normalizing constant がパラメーター  \theta の関数になっているためです。

Double Metropolis-Hastings Sampler

そこで Mele (2017) は、Liang (2010) で提案されている Double Metropolis-Hastings 法 を使い、パラメーターの事後分布を求めることを提案しました。 ここで double となっているのは、ネットワーク、パラメーター両方に対して Metorplis-Hastings 法を使うことに由来します。

まず、ネットワークをシミュレーションするアルゴリズムを書き下します。

Algorithm 1 - Metropolis-Hastings for Network Simulations

  1. パラメーター  \theta を固定する。イテレーション  r におけるネットワークが  g_{r} とする。
  2. 提案分布  q(\cdot| g_{r}) から からネットワーク  g^{\prime} をドローする。
  3. Acceptance ratio  \alpha(g_{r}, g^{\prime}) :=  \min\left( 1,  \frac{ \exp(Q(g^{\prime}, \theta)) }{ \exp(Q(g_{r}, \theta)) } \frac{ q(g_{r} | g^{\prime}) }{ q(g^{\prime} | g_{r}) } \right) を計算する。
  4. 一様分布  U[0, 1] から  u をドローする。
  5.  u \leq  \alpha(g_{r}, g^{\prime}) ならば、 g_{r} \leftarrow g{\prime} とする。
  6. 2に戻る

Acceptance ratio には normalizing constant が含まれていないので、計算が可能となります。 パラメーター  \theta のときの target distribution は


\begin{align}
\pi(g | \theta) =  \frac{\exp(Q(g, \theta))}{c(\mathcal{G} ,\theta)}
\end{align}

で、ネットワークをドローするにあたり  \theta を固定することで


\begin{align}
\frac{ \pi(g^{\prime} | \theta) }{ \pi(g_{r} | \theta) } \frac{ q(g_{r} | g^{\prime}) }{ q(g^{\prime} | g_{r}) }
= \frac{ \exp(Q(g^{\prime}, \theta)) }{ c(\mathcal{G} ,\theta) }  \frac{ c(\mathcal{G} ,\theta) } { \exp(Q(g_{r}, \theta)) } \frac{ q(g_{r} | g^{\prime}) }{ q(g^{\prime} | g_{r}) }
= \frac{ \exp(Q(g^{\prime}, \theta)) }{ \exp(Q(g_{r}, \theta)) } \frac{ q(g_{r} | g^{\prime}) }{ q(g^{\prime} | g_{r}) }
\end{align}

のように normalizing constant がキャンセルアウトされるためです。 これまでの連載でネットワークをシミュレーションする際に使っていたアルゴリズムは、この Algorithm 1 をです。 ネットワークの収束を改善するため、proposed network を生成する際に dyad の状態を一度に複数変更する large sampler を導入したのでした。

Double Metropolis-Hastings 法では、ネットワーク  g を補助変数として用いながら、パラメーター  \theta の事後分布を求めます。

Algorithm 2 - Double-Metropolis Hastings Sampler

  1. ネットワークのシミュレーション回数を  R とする。各イテレーション t におけるパラメーターを  \theta_{t} 、観察されたネットワークを  g とする。
  2. パラメーター  \theta の初期値  \theta_{0} を定める。
  3. 提案分布  q(\cdot| \theta) から  \theta^{\prime} をドローする。
  4. 観察されたネットワーク  g を初期値、パラメーターを  \theta^{\prime} として、Algorithm 1 を  R 回走らせる。 R 回後のシミュレーションにおけるネットワークを  g^{\prime} とする。
  5. Acceptance ratio 
\alpha(\theta, \theta^{\prime}, g^{\prime}, g)
:= \min\left(
1, \frac{ \exp(Q(g^{\prime}, \theta)) }{ \exp(Q(g, \theta)) }  \frac{ p(\theta^{\prime}) }{ p(\theta)  } \frac{ q(\theta | \theta^{\prime}) }{ q(\theta^{\prime} | \theta)} \frac{ \exp(Q(g, \theta^{\prime})) }{ \exp(Q(g^{\prime},  \theta^{\prime})) }
\right)
を計算する。
  6. 一様分布  U[0, 1] から  u をドローする。
  7.  u \leq  \alpha(\theta, \theta^{\prime}, g^{\prime}, g) ならば、 \theta \leftarrow \theta^{\prime} とする。
  8. 3に戻る。

以下、acceptance ratio がなぜこのような形で表せるかを確認します。

観察されたネットワークを  g とすると、今考える density は  \pi(\theta, \theta^{\prime}, g^{\prime}, | g) となります (ここでは補助変数として  g^{\prime} を考えていることに注意します)。 これは、

  1.  g のもとで  \theta を観察し ( p(\theta |g))
  2.  \theta を所与として  \theta^{\prime} を提案し ( q(\theta^{\prime} | \theta))
  3.  \theta のもとでネットワーク  g^{\prime} を観察する ( \pi(g^{\prime} | \theta^{\prime}))

と分解されますから、


\begin{align}
\pi(\theta, \theta^{\prime}, g^{\prime}, | g)
\propto p(\theta | g) q(\theta^{\prime} | \theta) \pi(g^{\prime} | \theta^{\prime})
= p(\theta) \pi(g, \theta) q(\theta^{\prime} | \theta) \pi(g^{\prime} | \theta^{\prime})
= p(\theta) \frac{  \exp(Q(g, \theta)) }{ c(\mathcal{G} ,\theta) } q(\theta^{\prime} | \theta) \frac{  \exp(Q(g^{\prime}, \theta^{\prime})) }{ c(\mathcal{G} ,\theta^{\prime}) }
\end{align}

と表せます。したがって 
\begin{align}
\frac{ \pi(\theta, \theta^{\prime}, g^{\prime}, | g) }{ \pi(\theta^{\prime}, \theta, g, | g^{\prime}) }
= \frac{ p(\theta) }{ p(\theta^{\prime})  } \frac{  \exp(Q(g, \theta)) }{ c(\mathcal{G}, \theta) }  \frac{  c(\mathcal{G}, \theta^{\prime}) }{ \exp(Q(g, \theta^{\prime})) }
\frac{  \exp(Q(g^{\prime}, \theta^{\prime})) }{ c(\mathcal{G} ,\theta^{\prime}) } \frac{ c(\mathcal{G} ,\theta) }{ \exp(Q(g, \theta)) }
\frac{q(\theta^{\prime} | \theta) }{q(\theta | \theta^{\prime}) }
\end{align}

\begin{align}
= \frac{ \exp(Q(g^{\prime}, \theta)) }{ \exp(Q(g, \theta)) }  \frac{ p(\theta^{\prime}) }{ p(\theta)  } \frac{ q(\theta | \theta^{\prime}) }{ q(\theta^{\prime} | \theta)} \frac{ \exp(Q(g, \theta^{\prime})) }{ \exp(Q(g^{\prime},  \theta^{\prime})) }
\end{align}

となり、acceptance ratio に normalizing constant が出てきません。

Algorithm 2 のポイントは、ネットワークのシミュレーションを行うときに、ネットワークの初期値を観察されたネットワークに設定していることです。 これは  \theta を固定したとき、ネットワークの事後分布が定常分布に収束することを保証します。詳しくは Mele (2017) の Appendix B を参照ください。

また、提案されたパラメーターの値によっては、ネットワークの収束に指数時間かかるケースもありえます。 そうした場合でも素早く収束するように large sampler を導入することが有効なのでした。

Double Metropolis-Hastings sampler を可視化すると以下のようになります。

f:id:komatsuna4747:20220314084419p:plain
Source: Park and Haran (2018)

Estimation in practice

Algorithm 1, 2 を実際に実装して、パラメーターの推定がうまくいくか試してみましょう。

ここでも、以下の GitHub repository のパッケージを利用します。

github.com

Network simulation

まずはネットワークを Algorithm 1 を用いてシミュレーションしてみます。今回は以下の効用関数を考えます。


\begin{align}
u_{i}(g, \theta) = \frac{\alpha}{2} \sum_{j=1}^{n} g_{ij} + \beta \sum_{j=1}^{n} \sum_{k \neq i, j}^{n} g_{ij}  g_{jk} + \frac{\gamma}{2} \sum_{j=1}^{n} \sum_{k \neq i, j}^{n} g_{ij}  g_{jk} g_{ki}
\end{align}

効用関数の各項はそれぞれ、直接のつながりから得られる benefit、人気のある友人から得られる benefit、共通の友人から得られる benefit と考えることができます。

このとき、ネットワーク形成ゲームは以下のようなポテンシャル関数によってポテンシャルゲームとして表現できます。


\begin{align}
Q(g, \theta) = \alpha \left(\frac{1}{2} \sum_{i=1}^{n} \sum_{j=1}^{n} g_{ij} \right) + \beta  \left( \frac{1}{2}\sum_{i=1}^{n} \sum_{j=1}^{n} \sum_{k \neq i, j}^{n} g_{ij}  g_{jk} \right) + \gamma  \left( \frac{1}{6}  \sum_{i=1}^{n} \sum_{j=1}^{n} \sum_{k \neq i, j}^{n} g_{ij}  g_{jk} g_{ki} \right)
\end{align}

ここで括弧でくくった項はそれぞれ、ネットワークの edge, twostar, triangle の数と一致します。

Algorithm 1 を計算するためには各イテレーションで

  1. ネットワーク  g' を提案分布からドローする
  2. accceptance ratio  \alpha(g_{r}, g^{\prime}) :=  \min\left( 1,  \frac{ \exp(Q(g^{\prime}, \theta)) }{ \exp(Q(g_{r}, \theta)) } \frac{ q(g_{r} | g^{\prime}) }{ q(g^{\prime} | g_{r}) } \right) を計算する

の2つが必要です。

まず今考える sampler では、  \frac{ q(g_{r} | g^{\prime}) }{ q(g^{\prime} | g_{r}) } = 1 となります。 ループのない任意のundirected network  g を考えます。そこからリンクを1つだけランダムに追加/削除するような sampler を考えると、


\begin{align}
q(g_{r} | g^{\prime}) = q(g^{\prime} | g_{r}) =  \frac{1}{N (N - 1) / 2}
\end{align}

となります。詳細は省きますが、large sampler を使う場合も   \frac{ q(g_{r} | g^{\prime}) }{ q(g^{\prime} | g_{r}) } = 1 が成立します。

そしてaccceptance ratio を計算していくわけですが、オーバーフローを回避するため


\begin{align}
\log(\alpha(g_{r}, g^{\prime})) =  \min\left(0,  Q(g^{\prime}, \theta)) - Q(g_{r}, \theta)) \right)
\end{align}

と自然対数をとります。ここで、ポテンシャル関数の各項がネットワークの統計量を表していたことを思い出せば、

Q(g^{\prime}, \theta)) - Q(g_{r}, \theta)) = \alpha * \Delta \text{edges} + \beta * \Delta \text{twostars} + \gamma * \Delta \text{triangles}

のように書けます。  \Delta \text{edges} \Delta \text{twostars} などのネットワークの統計量の変化分を change statistics と呼んだりします。

そして各イテレーションでU[0, 1] から乱数を発生させて

 \log(u) \leq \log(\alpha(g_{r}, g^{\prime}))

であれば  g_{r} \leftarrow g' のようにネットワークを更新していきます。

実際にネットワークをシミュレーションしてみます。

devtools::install_github("komatsuna4747/myergm")

library(myergm)
library(ggplot2)
set.seed(334)

# Number of nodes
N <- 200
# True parameters
theta_true <- c(-2.5, 1/N, 1/N)

# Simulate a fake network
net <-
  simulate_network(
    adjacency_matrix = matrix(0, nrow = N, ncol = N),
    coefEdges = theta_true[[1]],
    coefTwostars = theta_true[[2]],
    coefTriangle = theta_true[[3]],
    MCMC_interval = 1,
    MCMC_burnin = 0,
    MCMC_samplesize = 2[f:id:komatsuna4747:20220314084838p:plain]00000,
    p_one_node_swap = 0.001,
    p_large_step = 0.001,
    p_invert = 0.001,
    lambda = 0.2,
    level = "debug"
  )

今は空のネットワークから出発して、20万回イテレーションを回しました。ネットワークのシミュレーションが収束しているか確認してみます。

# Check convergence
df_stat_truenet <-
  tibble::tibble(
    iter = rep(1:nrow(net[[2]]), 3),
    stats = c(net[[2]][, 1], net[[2]][, 2], net[[2]][, 3]),
    terms = rep(c("edges","twostars", "triangles"), each = nrow(net[[2]]))
  )

g_stat_truenet <-
  df_stat_truenet %>%
  ggplot(aes(x = iter, y = stats, group = terms, color = terms)) +
  geom_line()

plot(g_stat_truenet)

f:id:komatsuna4747:20220314084838p:plain

空のネットワークから出発しても、10万回ぐらいイテレーションを回せばシミュレーションが収束してくれるようです。 したがって、後にパラメーターの事後分布を求めるときにネットワークをシミュレーションする回数を  R = 100,000 とおくことにします。

20万回イテレーションを回したあとに発生したネットワークを、今回は観察されたネットワークとして使用することにします。

Posterior distribution estimation

パラメーターの prior distribution を  \mathcal{N} (0, 10^2I) とおきます。 ここではランダムウォーク型 Metropolis-Hastings 法を使うことにし、パラメーター  \theta の提案分布として、 \mathcal{N} (\theta^{(t)}, \Sigma) を用いることにします。  \Sigma は最初は適当に  0.01 I とでもおいておきます。

まず、提案されたパラメーター  \theta' を用いて、観察されたネットワークを初期値としてネットワークを生成します。 その生成されたネットワークと観察されたネットワークから change statistics を計算し、acceptance ratio を計算します。

今考えている提案分布  q(\cdot| \theta) は対称ですから、  q(\theta' | \theta) / q(\theta | \theta') = 1 となります。したがって acceptance ratio は


\log \alpha(\theta, \theta^{\prime}, g^{\prime}, g)
:= \min\left(
0,  (\theta - \theta')^{\top} (\text{change stats}) - \frac{1}{2 * 10^2} \left( \Vert \theta' \Vert^2 - \Vert \theta \Vert^2 \right)
\right)

と書けます。

実際に Double Metropolis-Hastings 法を走らせてみましょう (推定は時間がかかります)。

# Prior
sigma0 <- diag(10, nrow = 3)
theta0 <- MASS::mvrnorm(n = 1, mu = c(0, 0, 0), Sigma = sigma0)
covpars <- diag(c(0.01, 0.01, 0.01), nrow = 3)

# Run DMH twice.
dmh_first <-
  run_DMH_cpp(
  adjacency_matrix = net[[1]],
  sigma_proposal = covpars,
  theta_init = theta0,
  sigma_prior = sigma0,
  terms_included = c(1, 1, 1),
  param_sample = 5000,
  MCMC_network_interval = 100000,
  level = "info"
)

結果を可視化してみます。

df_dmh_first <- 
  tibble::as_tibble(dmh_first)
colnames(df_dmh_first) <- c("edges", "twostar", "triangle")

plot_stat <- function(data) {
  data_tidy <-
    data %>%
    dplyr::mutate(iter = 1:nrow(.)) %>%
    tidyr::pivot_longer(-iter)
  
  g <-
    data_tidy %>%
    dplyr::mutate(name = forcats::fct_relevel(name, "edges", "twostar", "triangle")) %>%
    ggplot(aes(x = iter, y = value, group = name, color = name)) +
    geom_line() +
    facet_grid(. ~ name)
  
  return(g)
}

plot_stat(df_dmh_first)

f:id:komatsuna4747:20220314121459p:plain

うまく収束していないようです。提案分布の分散の設定がよくないのでしょう。 最初の1,000回分を捨てて、残りのデータから分散行列を計算し、それを使って Double Metropolis-Hastings をもう一度走らせてみます。

theta <-
  df_dmh_first %>%
  dplyr::slice(1000:nrow(.)) %>%
  tidyr::pivot_longer(dplyr::everything()) %>%
  dplyr::mutate(name = forcats::fct_relevel(name, "edges", "twostar", "triangle")) %>%
  dplyr::group_by(name) %>%
  dplyr::summarise(mean = mean(value)) %>%
  dplyr::ungroup() %>%
  dplyr::pull(mean) %>%
  as.matrix()
  
covpars <- cov(df_dmh_first %>% dplyr::slice(1000:nrow(.)))

dmh_second <-
  run_DMH_cpp(
    adjacency_matrix = net[[1]],
    sigma_proposal = covpars,
    theta_init = theta,
    sigma_prior = sigma0,
    terms_included = c(1, 1, 1),
    param_sample = 10000,
    MCMC_network_interval = 100000,
    level = "info"
  )

推定結果を確認します。

df_dmh_second <- 
  tibble::as_tibble(dmh_second)
colnames(df_dmh_second) <- c("edges", "twostar", "triangle")

plot_stat(df_dmh_second)

f:id:komatsuna4747:20220314121959p:plain

うまく収束しているように見えます。

最初の1,000回分を burnin として捨て、残りのデータから要約統計量を計算します。

summary(df_dmh_second %>% dplyr::slice(1000:nrow(.)))
     edges           twostar             triangle        
 Min.   :-3.285   Min.   :-0.008715   Min.   :-0.161898  
 1st Qu.:-2.798   1st Qu.: 0.006354   1st Qu.:-0.050359  
 Median :-2.661   Median : 0.011002   Median :-0.027502  
 Mean   :-2.657   Mean   : 0.010702   Mean   :-0.030452  
 3rd Qu.:-2.511   3rd Qu.: 0.015185   3rd Qu.:-0.007288  
 Max.   :-2.017   Max.   : 0.029968   Max.   : 0.045579  

真のパラメーターを  (\text{edges}, \text{twostar}, \text{triangle}) = (-2.5, 0.005, 0.005) と設定していたので、triangle に関しては推定がうまくいってないように見えます。 triangle の係数の推定は難しいと言われるようですが、また別途機会を見つけて原因を解明したいです。

Concluding remarks

今回は Mele (2017) の提案に沿って、Double Metropolis-Hastings 法により ERGM のパラメーターの推定を試みました。

ERGM の推定でよく使われる R のパッケージは {ergm} ですが、そのパッケージでは Double Metropolis-Hastings ではなく、Monte Carlo for Maximum Likelihood Estimation が使われています。 これはネットワークの Markov chain を生成して、そこから複数のネットワークをサンプリングし ERGM のパラメーターを推定するというものです。 詳細は Statnet チームのドキュメント、または弊社研究員前嶋のブログが参考になるでしょう。

statnet.org

meana0.hatenablog.com

これまでは ERGM の理論的な話が多かったので、次回は ERGM の応用例を見ていくことにします。

References

  • Liang, F. (2010). A double Metropolis–Hastings sampler for spatial models with intractable normalizing constants. Journal of Statistical Computation and Simulation, 80(9), 1007-1022.
  • Mele, A. (2017). A structural model of dense network formation. Econometrica, 85(3), 825-850.
  • Park, J., & Haran, M. (2018). Bayesian inference in the presence of intractable normalizing functions. Journal of the American Statistical Association, 113(523), 1372-1390.
  • 鎌谷研吾. (2020). モンテカルロ統計計算. 講談社.


▼本連載のほかの記事はこちら buildersbox.corp-sansan.com

© Sansan, Inc.