やったこと
テンポラル・ネットワーク（構造が時間とともに変化するネットワーク）の分析モデルの一つ temporal fitness model (Kobayashi et al., 2019) に登場する以下の非線形連立方程式をPythonで実装しました。
\tag{8}H_i(\boldsymbol{a^*})\equiv\sum_{j:j\neq i}\frac{m_{ij}^o-\tau a_i^*a_j^*}{1-a_i^*a_j^*}=0,~~~~~~~ \forall\,i=1,...,N
(8)式の簡単な説明
Temporal fitness model ではネットワークを構成する各ノード $i\in\{1,\dots,N\}$ は、活動量（activity level）$a_i\in(0,1]$という固有の値を持っていると仮定します。すると、一定の時間の間にノード$i$とノード$j$がつながる確率$u(a_i, a_j)$を、それぞれの活動量の積$a_ia_j$で定義できます。これは2国間の貿易量がそれぞれの国の経済規模に比例する重力モデルに似た定式化です。
上で定義された接続確率を用いると、あるノードのペアが一定時間に何回接続するかを二項分布に基づいて計算することができ、その理論的な分布を使って統計的検定が可能になります。ただ、そのためにはノードの活動量をデータから推定する必要があり、その過程で式(8)の非線形連立方程式が登場するわけです。式の中の$a_i^*$は最尤推定値で$m_{ij}^o,~\tau$は外生変数（パラメータ）です。詳細は (Kobayashi et al., 2019) を参照してください。
Pythonで(8)式をどう表現するか
式の表現には、NumPyを使いました。ポイントは以下の通りです：
- ノード$i$と$j$の活動量の積$a_ia_j$は行列（列ベクトル@行ベクトル）で表現できる
- 個々の非線形な項は、NumPy配列の２種類の積（*と@）を使えば表現できる
- $i=1,…N$の一本一本の式は
numpy.ndarray.sum(axis=1)（行和）で表現できる - 各項に出てくる$m_{ij}^o$は、対角成分が全てゼロの行列の(i,j)成分
- 添え字が一致する項は和に入れない
実装
ソルバーにはSciPyのoptimizeモジュールのrootを使いました。詳しい利用方法はドキュメントを参照してください。ソルバーを使って解く関数を、上記のポイントにしたがって次のように定義します。aamは$m_{ij}^o$を要素にもつ$N \times N$の行列（NumPyの2次元配列）です。
def h(a):
a = np.array(a)
monomials = (aam - tau*(a.reshape(-1,1)@a.reshape(1,-1))) / (np.ones((n,n)) - (a.reshape(-1,1)@a.reshape(1,-1)))
diag_ = np.diag_indices(n)
monomials[diag_] = 0
poly = monomials.sum(axis=1)
return poly
このように定義したh(a)を初期値とともにscipy.optimize.rootに渡せば解くことができます。
a = optimize.root(h, a_init, method = "krylov")
おわりに
ここまで読んでいただきありがとうございます。もともとこの連立方程式の求解問題は、卒業研究に取り組んでいた時に直面したものでした。Temporal fitness model を使って何をしたのかに興味がある方がいらっしゃれば、こちらをご覧いただけると嬉しいです。
動的ネットワーク分析を用いたInstagram上のトレンド推定に関する研究 / Instagram Trend Analysis - Speaker Deck
また、Temporal fitness model (ST filter)自体を使ってみたい方は、ゼミの先輩が先月パッケージST_filterをリリースされたので、ぜひ活用してください。
参考文献
- Kobayashi, Teruyoshi, Taro Takaguchi, and Alain Barrat. 2019. “The structured backbone of temporal social ties.” Nature communications 10 (1): 1–11.
- 非線形連立方程式ってPythonで簡単に解けるんですね - Qiita
- 経済学のためのPython入門