f:id:shiori_baba:20220316114139j:plain

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

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

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 法の手順は以下の通りです。

パラメーター $\theta$ の初期値 $\theta_{0}$ を定める。
提案分布 $q(\cdot| \theta)$ から $\theta^{\prime}$ をドローする。
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)$ を計算する。
一様分布 $U[0, 1$ ] から $u$ をドローする。
$u \leq \alpha(\theta, \theta^{\prime})$ ならば、 $\theta \leftarrow \theta^{\prime}$ とする。
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

パラメーター $\theta$ を固定する。イテレーション $r$ におけるネットワークが $g_{r}$ とする。
提案分布 $q(\cdot| g_{r})$ からからネットワーク $g^{\prime}$ をドローする。
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)$ を計算する。
一様分布 $U[0, 1$ ] から $u$ をドローする。
$u \leq \alpha(g_{r}, g^{\prime})$ ならば、 $g_{r} \leftarrow g{\prime}$ とする。
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

ネットワークのシミュレーション回数を $R$ とする。各イテレーション $t$ におけるパラメーターを $\theta_{t}$ 、観察されたネットワークを $g$ とする。
パラメーター $\theta$ の初期値 $\theta_{0}$ を定める。
提案分布 $q(\cdot| \theta)$ から $\theta^{\prime}$ をドローする。
観察されたネットワーク $g$ を初期値、パラメーターを $\theta^{\prime}$ として、Algorithm 1 を $R$ 回走らせる。 $R$ 回後のシミュレーションにおけるネットワークを $g^{\prime}$ とする。
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)$ を計算する。
一様分布 $U[0, 1$ ] から $u$ をドローする。
$u \leq \alpha(\theta, \theta^{\prime}, g^{\prime}, g)$ ならば、 $\theta \leftarrow \theta^{\prime}$ とする。
3に戻る。

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

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

$g$ のもとで $\theta$ を観察し ( $p(\theta |g)$ )
$\theta$ を所与として $\theta^{\prime}$ を提案し ( $q(\theta^{\prime} | \theta)$ )
$\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 を計算するためには各イテレーションで

ネットワーク $g'$ を提案分布からドローする
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 Tech Blog

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

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