ネットワークグラフの中心性

前の記事でサンプルとして(折れ線グラフや棒グラフではなく、グラフ理論で言うところの)グラフを扱ったのでもう少し紹介します。
参考: Matplotlibの配色を別の処理でも流用したい

グラフについて分析するとき「どの点が重要なのか」という問いは非常に自然なものです。
そして、「中心にある点が重要なんじゃないか」と考えることもそこそこ自然な発想になります。
ただ、扱う対象が普通の幾何学的な図形ではなく、点とそれらの間にあるつながりという抽象的な概念で構成されたグラフの場合、
どの点が中心なのかというのは非自明な問題になります。
可視化してみれば確かにどれかの点が真ん中らへんにあるように見えるのですが、それは単に可視化の際にたままたそこにプロットされたというだけで、
グラフ理論で言うところの中心ではないからです。

この問題に対応するために、複数の中心が定義され、中心らしさを表す指標がいくつか提案されています。
networkxにも色々定義されているので、その中でも定義が単純でわかりやすいものを紹介します。
networkxのドキュメントではこちらに該当します。
Centrality

今回の記事では次のグラフをサンプルに使います。
何か特別な名前がついてるやつではなく、僕が今適当に構成したものです。


import matplotlib.pyplot as plt
import networkx as nx

G = nx.Graph()
G.add_edges_from([
    (0, 1),
    (1, 2),
    (2, 3),
    (3, 0),
    (1, 3),
    (0, 4),
    (4, 5),
    (5, 0),
    (5, 6),
])

# 可視化
fig = plt.figure(facecolor="w")
ax=fig.add_subplot(1, 1, 1, title="Sample Graph")
nx.draw_networkx(
    G,
    node_color="c"
)

さて、サンプルが用意できたので、次数中心性/近接中心性/媒介中心性の3種類の中心性の定義を紹介します。

次数中心性
一番単純でわかりやすいのが次数中心性(degree centrality)です。
これは、各ノードの次数(そのノードにつながっている、辺の本数)を指標とするものです。
正確には、各ノードの次数は最大でも「ノード数-1」までの値しか取れないので、「ノード数-1」で割って正規化したものを使います。
networkxではdegree_centralityで計算できます。


print("ノード番号, 次数中心性")
for k, v in nx.degree_centrality(G).items():
    print(k, ",", v)
"""
ノード番号, 次数中心性
0 , 0.6666666666666666
1 , 0.5
2 , 0.3333333333333333
3 , 0.5
4 , 0.3333333333333333
5 , 0.5
6 , 0.16666666666666666
"""

計算してみた各値を6(=7-1)倍すると、全部整数になり、各ノードから伸びてる辺の数に一致することがわかります。

近接中心性
距離がベースになっていて、円の中心などとイメージが近いのが近接中心性(closeness centrality)です。
これは、ある点から、そのほかの全ての点へのそれぞれの距離の平均値を元に決まります。
ただ、ほかの全ての点に近い(つまり距離が小さい)ほど、中心と見なしたいので、
距離の平均の逆数を取ります。
今回のグラフでは各辺に重み付けなどしていない(つまり全部重さ1とみなす)ので、ここで言う距離というのは最短経路を通った時に通過する辺の数です。
networkxではcloseness_centralityで計算できます。

※ 2020/12/24 コードに誤りがあることを教えていただき修正しました。


print("ノード番号, 近接中心性")
for k, v in nx.closeness_centrality(G).items():
    print(k, ",", v)
"""
ノード番号, 近接中心性
0 , 0.75
1 , 0.6
2 , 0.42857142857142855
3 , 0.6
4 , 0.5454545454545454
5 , 0.6
6 , 0.4
"""

試しに一つ具体的に計算しておきましょう。
0番のノードに着目します。
このノードの1~6番目までのノードの距離はそれぞれ、$1,2,1,1,1,2$ で、合計は$8$です。なので平均は$4/3$になります。
これの逆数をとったものが$3/4=0.75$なので、ライブラリの結果と一致しました。

媒介中心性
最後が媒介中心性(betweenness centrality)です。
個人的には、この記事で紹介した3種類の中では一番よく使います。
(主観的な感想です。媒介中心性がほかの二種に比べて理論的に優れているということを主張するものではありません。)
これは、着目している点以外の2点を結ぶ最短経路のうち、その点を通過するものの割合です。
networkxではbetweenness_centralityで計算できます。


print("ノード番号, 媒介中心性")
for k, v in nx.betweenness_centrality(G).items():
    print(k, v)
"""
ノード番号, 媒介中心性
0 0.6
1 0.13333333333333333
2 0.0
3 0.13333333333333333
4 0.0
5 0.3333333333333333
6 0.0
"""

さて、これは個別に解説したいと思います。
まず、2, 4, 6 の媒介中心性は0になっています。
これは、可視化した図を見るとわかるのですが、それぞれの点について、ほかの2点をどう選んでも、対象の点を通ると遠回りになるからです。

次にわかりやすいのは
5番です。まず、5番以外の6点から2点を選ぶ方法は全部で$6*5/2=15$通りあります。
その中で、5番を経由すると最短になるのは、6番の点と、そのほかの5点のどれかを結ぶ5通りです。
そのため、5版のノードの媒介中心性は$5/15=1/3$になります。
0番も同様に数えると、4,5,6のどれかと、1,2,3のどれかを結ぶ$3*3=9$通りの最短経路が0番を通過し、
媒介中心性は$9/15=3/5=0.6$になります。

さて、のこるは同じ値になっている1番と3番です。
まずは1番の方ですが、最短経路で1番のノードを通るのは2番のノードと、0,4,5,6のどれかを結ぶ4本になります。
となると、媒介中心性は$4/15=2.666…$になっても良さそうなのですが、実際にはこの$1/2$の値になっています。
実はこれは2番のノードと、0,4,5,6のどれかを結ぶ最短経路がそれぞれ複数あるのが原因です。
1番のノードの代わりに3番のノードを通っても良いわけです。
ということで、1番を通る最短経路を数える時に$0.5$倍して数えるので上の計算例の値になります。
3番も同様です。

辺の媒介中心性
さて、ノードの中心性を3つ紹介しましたが、最後の媒介中心性はノードだけではなく、辺についても定義できます。
ノードの場合は自分以外の2点でしたが、辺の場合はグラフ中の全てのノードから2点選び、
それらの最短経路のうち何本が自身を通過するかで定義されます。
networkxではedge_betweenness_centralityで計算できます。

計算結果が少数が複雑で、ぱっと見わかりにくかったので、
ノードの組みの場合の数、つまり$7*6/2=21$通りを掛けて、整数化したものを表示してみました。


print("辺, 辺の媒介中心性*21")
for k, v in nx.edge_betweenness_centrality(G).items():
    print(k, v*21)
"""
辺, 辺の媒介中心性*21
(0, 1) 6.0
(0, 3) 6.0
(0, 4) 4.0
(0, 5) 8.0
(1, 2) 3.0
(1, 3) 1.0
(2, 3) 3.0
(4, 5) 2.0
(5, 6) 6.0
"""

(1, 3) が 1 なのは、まさにその1と3をつなぐ時以外は通らないからですね。
一方で(5, 6)は、6とほかの点をつなぐ時は必ず通らないといけないので、中心性が高くなっています。

可視化した図では(5, 6)のエッジはかなり端っこにあるように見えるのですが、
媒介中心性は2位タイの高さということで、
可視化した時に真ん中にあるかどうかと、グラフ理論で言うところの中心は違うと言うことの雰囲気が伝わればと思います。

最後に、もう少し大きめのグラフで、それぞの中心性を算出して比較してみましょう。
計算した中心性はノードのサイズで示しました。


# ランダムにグラフデータ生成
while True:
    G=nx.random_graphs.fast_gnp_random_graph(25, 0.1)
    if nx.is_connected(G):
        # 連結なグラフができたら抜ける
        break


# 3種類の中心性をそれぞれ計算
dc_dict = nx.degree_centrality(G)
dc_list = np.array([dc_dict[n] for n in G.nodes])

cc_dict = nx.closeness_centrality(G)
cc_list = np.array([cc_dict[n] for n in G.nodes])

bc_dict = nx.betweenness_centrality(G)
bc_list = np.array([bc_dict[n] for n in G.nodes])

# 中心性をノードのサイズで表して可視化
pos = nx.spring_layout(G, iterations=300)
fig = plt.figure(facecolor="w", figsize=(8, 24))

ax=fig.add_subplot(3, 1, 1, title="Degree centrality ")
nx.draw_networkx(
    G,
    pos=pos,
    ax=ax,
    node_color="c",
    node_size=dc_list*300/np.mean(dc_list),
)
ax=fig.add_subplot(3, 1, 2, title="Closeness centrality")
nx.draw_networkx(
    G,
    pos=pos,
    ax=ax,
    node_color="c",
    node_size=cc_list*300/np.mean(cc_list),
)
ax=fig.add_subplot(3, 1, 3, title="Betweenness centrality")
nx.draw_networkx(
    G,
    pos=pos,
    ax=ax,
    node_color="c",
    node_size=bc_list*300/np.mean(bc_list),
)
plt.show()

出力された結果がこちら。

Matplotlibの配色を別の処理でも流用したい

Python(と jupyter notebook)でデータを可視化する場合、色を16進法のRGBで指定できるライブラリは多くあります。
Matplotlibがベースになっているものは、そのカラーマップを指定できることも多いのですし、
「rは赤」、「bは青」など一部の色はアルファベットや色名で指定できるのですが、
もっと多くの色を使いたかったり、値によってグラデーションをつけたい場合で逐一RGBを構築するのは結構な手間です。

そこで、Matplotlibの配色をそのまま流用できないかと思って調べてみました。
結論から言うと、結構簡単に使えそうです。

まず、配色そのもののデータは、
matplotlib.cm と言うモジュールに含まれています。
配色はその名前で指定しますが、名前と実際の色の対応はこちらのリファレンスをみると良いでしょう。
Colormap reference

使いたいカラーマップが決まったら、cm.get_cmap() か、 cmの属性として、使うことができます。
要するに次の2行は同じものです。


cm.get_cmap("Greens")
cm.Greens

さて、どちらもカラーマップのオブジェクトを返してくれますが、
そのカラーマップのオブジェクトにに数値を渡すと、RGBのタプルを返してくれます。


import matplotlib.cm as cm

print(cm.get_cmap("Greens")(0.7, alpha=0.5))
# (0.18246828143021915, 0.5933256439830834, 0.3067589388696655, 0.5)
print(cm.get_cmap("Paired")(3))
# (0.2, 0.6274509803921569, 0.17254901960784313, 1.0)

渡す数値ですが、連続的に色が変化するものには、 0〜1の値を渡します。
色の値が不連続な(要はリファレンスで、Qualitativeのカテゴリにあるもの)は、0〜1の値で渡しても大丈夫ですが、
整数値で0,1,2,3などを指定しても大丈夫です。
これらは、1と1.0や2と2.0など、同じ値でも整数型と浮動小数型で結果が変わるので注意してください。
ちなみに、値はリスト形式で複数同時に渡しても大丈夫です。

さて、最初の話に戻りますが、このRGB値のタプルを他のライブラリ等で使うには、16進法の文字列に変換する必要があります。
255倍して16進法の文字列に変化して、シャープをつけて結合するコードを自分で書いてもいいのですが、
なんと Matplotlibにその関数が用意されていました。

matplotlib.colors.rgb2hex です。
これはなぜか、色のリストは受け取ってくれないので、順番に適用していかないといけないのですが、
RGBのタプルを16進法文字列に手軽に変換してくれます。
(keep_alpha=Trueを指定すると透明度も含めてくれます。デフォルトはFalseです。)

試しにカラーマップから10色取り出してみましょう。


import matplotlib.colors as mcolors
import matplotlib.cm as cm
import numpy as np

cmap = cm.get_cmap("BuGn")
for rgb in cmap(np.arange(0, 1, 0.1)):
    print(mcolors.rgb2hex(rgb))

"""
#f7fcfd
#e9f7fa
#d6f0ee
#b8e4db
#8fd4c2
#65c2a3
#48b27f
#2f9858
#157f3b
#006428
"""

10個の色が取り出せましたね。

最後に何か例を出しておきたいので、networkx で作成したグラフに中心性で色をつけてみました。
(グラフの中心性には複数の種類がありますが、今回は媒介中心性を使いました。)


import networkx as nx
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as mcolors


while True:
    # ランダムにグラフを生成する
    G = nx.random_graphs.fast_gnp_random_graph(15, 0.2)
    # 連結なグラフが生成できたらループを抜ける
    if nx.is_connected(G):
        break

# 媒介中心性の計算
centrality = nx.betweenness_centrality(G)
# 辞書形式なので、ノードの順番と揃えてリスト化する。
centrality_list = np.array([centrality[node] for node in G.nodes])
# 媒介中心性を0〜1に正規化する
color_level = centrality_list - min(centrality_list)
color_level/=max(color_level)
# ノードの色の生成
rgb_list = cm.get_cmap("Oranges")(color_level, alpha=0.8)
node_color=[mcolors.rgb2hex(rgb, keep_alpha=True) for rgb in rgb_list]

# グラフの可視化
fig = plt.figure(figsize=(8, 8), facecolor="w")
ax = fig.add_subplot(1,1,1)
nx.draw_networkx(
                G,
                node_color=node_color,
                node_size=500,
                edge_color="#aaaaaa",
                node_shape="s"
            )

出力がこちら。

媒介中心性が高いところが色が濃くなっているのがわかります。

WordPressにお問い合わせフォームを作る

このブログに問い合わせフォームをつけることにしたのでその手順のメモです。
複数のプラグインがあるようですが、比較的情報が多かった Contact Form 7 を使うことにしました。

まずはプラグインをインストールします。
– ダッシュボードの左ペインのプラグインをクリック
– Contact Form 7 を検索
– 今すぐインストールボタンからインストール
– インストールしたら有効化ボタンをクリック
– 左ペインに お問い合わせ メニューができたことを確認

新しくできた お問い合わせメニューを選択すると、最初からコンタクトフォーム1というフォームができていました。
設定項目がたくさんあるのですが、ざっと見た限りではそのまま使って問題なさそうなのでデフォルトで使います。

このショートコードをコピーして、投稿、固定ページ、またはテキストウィジェットの内容にペーストしてください:
というコメントの下に貼り付け用のコードがあるのでそれをコピーします。

これをブログ内のどこかに貼れば完成するのですが、サイドバーに置くことにしました。
ウィジェットの一覧から「テキスト」を選択し、ブログサイドバーに追加します。
そして、テキストの本文部分に先ほどの貼り付け用のコードを挿入して完成です。

あとは各ページに表示されたので送信テストを行いました。

DataFrameの日付の欠損行を埋める方法

日付に限らず連番等でも使える方法ですが、自分が日単位の時系列データで行うことが多いのでそれで説明します。

DBなどからデータを日単位で集計してpandasのDataFrameを作った時、集計対象のデータがなかった日は行ごと欠損してしまった状態になります。

例えば次のような感じです。


print(df.head(5))
"""
         date  value
0  2020-01-02      9
1  2020-01-06      3
2  2020-01-07      4
3  2020-01-09      9
4  2020-01-11      2
"""
print(df.tail(5))
"""
          date  value
15  2020-02-01      3
16  2020-02-02      2
17  2020-02-09      3
18  2020-02-11      2
19  2020-02-18      2
"""

このままで困らないこともあるのですが、累積和をとるときや、matplotlibで可視化するときなど、
欠損してる日付を補完しておきたいことがあります。

これまで、補完対象のDataFrameを別途構成してappendすることが多かったのですが、
必要な日付の一覧を持ったDataFrameと結合(SQLでいうJoin,pandasの関数ではMerge)すると手軽に補完できることに気づきました。

具体的には次のようなコードになります。


# dates に必要な期間の日付の一覧が入ってるとします。
date_df = pd.DataFrame({"date": dates})

# date_df と 結合する
df = pd.merge(date_df, df, on="date", how="left")
# NaNを 0で埋める
df[["value"]] = df[["value"]].fillna(0)

この例だと単純なので、不足している分のデータFrameを作ってたすのと比べて、あまりメリットを感じないのですが、
これが、例えば複数のキーに対してそれぞれ日付データを全部揃えたいケースになると、かなり楽になります。

例えば、元のデータフレームが次だったとします。


print(df2)
"""
    key        date  value
0  key1  2020-01-03      5
1  key1  2020-01-04      3
2  key1  2020-01-10      5
3  key2  2020-01-02      4
4  key2  2020-01-03      1
5  key2  2020-01-04      9
6  key3  2020-01-04      2
7  key3  2020-01-06      5
8  key3  2020-01-09      8
"""

このDataFrameの key1,key2,key3 に対して、 2020-01-01〜2020-01-10の行を全て揃えたいとします。
このようなときは、次鵜のようにして、keyの値と日付の値のペア全てのDataFrameを作ってそれと結合すると簡単に保管できます。


# key と 日付のペアを網羅したDataFrameを作る
keys, dates = np.meshgrid(
        ["key1", "key2", "key3"],
        [(datetime(2020, 1, 1) + timedelta(days=i)).strftime("%Y-%m-%d") for i in range(10)]
    )

key_date_df = pd.DataFrame(
        {
            "key": keys.ravel(),
            "date": dates.ravel(),
        }
    )

# 結合してソート
df2 = pd.merge(key_date_df, df2, how="left").sort_values(["key", "date"])
# NaNを 0で埋める
df2[["value"]] = df2[["value"]].fillna(0)
# インデックのリセット
df2.reset_index(inplace=True, drop=True)

途中で meshgridを使いましたが、meshgridに慣れてない場合は別の方法でも大丈夫です。

Pythonで連続した日付のリストを作る

日付の連番を文字列で必要になったので、Pythonで生成する方法を二つメモしておきます。

一つ目は、 標準ライブラリである datetime を使うものです。
開始日を生成して、必要な日数だけtimedeltaで差分を加算したものをリスト化したら得られます。
生成したリストはdatetime.datetime型なので、strftimeで文字列に変換して完成です。


from datetime import datetime, timedelta

# 日付のリスト生成()
date_list = [datetime(2020, 1, 25) + timedelta(days=i) for i in range(10)]
# 文字列に変換
date_str_list = [d.strftime("%Y-%m-%d") for d in date_list]
print(date_str_list)
"""
['2020-01-25', '2020-01-26', '2020-01-27', '2020-01-28',
'2020-01-29', '2020-01-30', '2020-01-31',
'2020-02-01', '2020-02-02', '2020-02-03']
"""

もう一つはpandasのdate_range関数を使います。
いくつかみて回った限りではこちらの方が人気のようです。
生成されるのが、DatetimeIndex なので、DataFrameのIndexで使いたい場合はこちらの方が便利なのだと思います。
また、生成するデータの頻度を指定するオプションが異常なほど充実しています。
参考: Time series / date functionality

とりあえず、同じデータを生成してみます。


import pandas as pd 

date_index = pd.date_range("2020-01-25", periods=10, freq="D")
print(date_index)
"""
DatetimeIndex(['2020-01-25', '2020-01-26', '2020-01-27', '2020-01-28',
               '2020-01-29', '2020-01-30', '2020-01-31', '2020-02-01',
               '2020-02-02', '2020-02-03'],
              dtype='datetime64[ns]', freq='D')
"""

# 配列に変換して必要な文字列に加工
date_ary = date_index.to_series().dt.strftime("%Y-%m-%d")
print(date_ary.values)
"""
['2020-01-25' '2020-01-26' '2020-01-27' '2020-01-28' '2020-01-29'
 '2020-01-30' '2020-01-31' '2020-02-01' '2020-02-02' '2020-02-03']
"""

これだけだと、ちょっと手間が余計にかかっていて、2つ目の方法にメリットがないように見えますが、
date_rangeは指定できる引数の種類が多く、場合によってはかなり柔軟に対応できます。

たとえば、開始日時と件数の代わりに、開始日時と終了日時で指定したり、終了日時とデータ件数で指定できます。
次の3行は全て同じ結果を返します。


pd.date_range("2020-01-25", periods=10, freq="D")
pd.date_range(start="2020-01-25", end="2020-02-03", freq="D")
pd.date_range(end="2020-02-03", periods=10,  freq="D")

また、時間単位や月単位、月単位といった頻度もfreqで指定できますが、
平日のみとか、毎月の15日と月末日など、datetimeで実装するには少し面倒なものも手軽に作れます。
再掲ですが、こちらのリファレンスを見ると色々あって面白いです。
Time series / date functionality

PyPIのパッケージをcondaでインストールする方法

前回の記事で、condaの公式リポジトリとconda-forgeを探して、それでも見つからなかったパッケージは環境を分けてpipで入れるという話を書きましたが、
一応、PyPIのパッケージをcondaで入れる方法は存在します。

それが、自分でcondaバッケージをビルドする方法で、(スムーズに進めば)具体的には次の手順でできます。
(前回の記事の反省を踏まえてさっさとコマンドを載せておきます。)


conda skeleton pypi [パッケージ名]
conda build [パッケージ名]
conda install --use-local [パッケージ名]

ドキュメントはこちら:Conda-build documentation

試しにこのブログを書くのに欠かせない pycodestyle_magic を入れてみます。これはcondaの公式リポジトリにもconda-forgeにも今日時点では存在しません。
準備として、必要なパッケージである、次の3つをcondaで入れておきます。
これがconda install の場合との大きな違いで、依存するライブラリは事前自分で入れておかないとエラーになります。
flit はこれが何かよく理解していないのですが、入れておかないとビルドする時に、 No module named ‘flit’ というエラーが出ました。


conda install pycodestyle
conda install flake8
conda install flit

さて、準備できたら次のコマンドを順番に実行します。


conda skeleton pypi pycodestyle_magic
conda build pycodestyle_magic
conda install --use-local pycodestyle_magic

conda skelton を実行した段階で、カレントディレクトリ直下に、meta.yamlファイルを含むディレクトリが出来上がります。


$ ls
pycodestyle_magic
$ ls pycodestyle_magic/
meta.yaml

次のconda build も成功するはずだったのですが、結局エラーが出たので対応します。
(3行目の conda installはまだ実行できず。)

出たエラーはこれ。


ModuleNotFoundError: No module named 'flit'

事前にインストールしているのに見つけられていないようです。
meta.yamlファイルを開いて、追記してみます。
(numpyで似たようなことをして解決している人がいたので真似しました。結果的に成功したようです。)


vim pycodestyle_magic/meta.yaml

#requirements の host: と run: に - flit を追記して保存。

requirements:
  host:
    - pip
    - python
    - flit
  run:
    - python
    - flit

さて、もう一回、conda build pycodestyle_magicすると成功しました。


####################################################################################
Resource usage summary:

Total time: 0:00:23.6
CPU usage: sys=0:00:00.3, user=0:00:00.6
Maximum memory usage observed: 75.1M
Total disk usage observed (not including envs): 264B


####################################################################################
Source and build intermediates have been left in /Users/[ユーザー名]/.pyenv/versions/[インストールしたAnacondaのディレクトリ]/conda-bld.
There are currently 1 accumulated.
To remove them, you can run the ```conda build purge``` command

メッセージにある通り、Anacondaのインストールフォルダ配下にビルドのアウトプットは移動しており、
さっきのmeta.yamlファイルがあった所の中身はそのままです。

conda-bld/ の下には、 pycodestyle_magic_1582538942427 というディレクトリができていて、よくわからないファイルがあります。

さて、buildが成功したので最後にインストールです。
(紛らわしい書き方して済みませんが、Anacondaのディレクトリではなく、skeletonやbuildしたのと同じディレクトリでそのまま)


# 今度は成功する。
conda install --use-local pycodestyle_magic

# インストール結果の確認
$ conda list | grep pycodestyle
pycodestyle               2.5.0                    py37_0
pycodestyle_magic         0.5                      py37_0    local

結果確認は、わかりやすいように普通にcondaで入れたpycodestyleも表示しましたが、
このskeletonを使う方法で入れたパッケージは local マークがつきます。

最後に後片付けです。
カレントディfレクトリに、skeletonでできた meta.yamlファイルが入ったディレクトリがあると思いますが、
インストールが終わったら不要なので消してよいと思います。
また、
conda build purge コマンドで、
[Anacondaのインストールディレクトリ]/conda-bld
にできていた一次ファイルたちを消せます。


rm -r [パッケージ名]
conda build purge

あとは入れたパッケージを使うだけです。
今回はjupyterのmagicコマンドだったので、
%load_ext pycodestyle_magic

%%pycodestyle
をそれぞれ試しました。

これで、conda-forgeにもないパッケージもPyPIからインストールできました。

condaの基本的な使い方

PyCon JP 2019 に参加して、
Anaconda環境運用TIPS 〜Anacondaの環境構築について知る・質問に答えられるようになる〜
という講演を聞いてから、自分の環境がcondaとpipの混在状態で結構リスキーだったことが懸念でした。
参考: 資料のページ

(実際何度も環境壊れていますし。)
会社で使っている端末では環境を極力condaのみで作り直していたのですが、実は自宅の端末では相変わらず混在していました。
そんな時に先日、SSDが壊れてMBPを修理に出す機会があり、完全に初期化されたので今回はできる限りpipを使わないように作り直しました。

まず、Anaconda入れるところまでは昔の記事の通りです。
参考: Macにpyenvとanacondaをインストールする

これまではこの後pip install [ライブラリ名]で、どんどんライブラリを入れていましたが、これはとても危険な行為です。
condaとpipの仕組みには互換性がないからです。

condaのドキュメントにも Using pip in an environment というセクションで注意事項が書かれています。
要約すると、pipを使うなら環境を分けることと、condaでできるだけ多くのライブラリをインストールした後に残りをpipで入れるべきということです。

免責事項
さて、ここからcondaの説明など書いていきますが、僕自身がパッケージの開発者でもなく、環境構築等を専門にするエンジニアでもないただのユーザーなので、
とりあえず自分はこういう理解で使っています、というレベルの内容になります。
厳密には各自公式のドキュメントを確認の上、自己責任でご利用お願いします。
(どちらかというと、他の記事の pip install 〜と書いてるところにこそ、この免責事項が必要ですね。)

さて、 pip と conda ですが、どちらもパッケージマネージャー(パッケージを管理するツール)です。
そして管理の方法が違い、互換性がないので混ぜると二重管理になって競合し、環境の破壊等が発生しえます。

インターネット上のリポジトリで管理されているパッケージを手元の端末に追加したりアップデートしたりできるのですが、
その大元になるリポジトリが違います。
pip は PyPI 、で、 condaはrepo.anaconda.comです。

condaの最大のデメリットが、PyPIに比べてrepo.anaconda.comで配布されているパッケージの少なさだと思います。
この点を補うために、condaでは別の場所(チャンネル)からもパッケージを追加できます。
特に conda-forge は比較的充実しているのでオススメです。
ドキュメントの中でも紹介されているので安心して使えるかと思います。
参考: Conda channels

余談ですが、Python初心者の頃、-c conda-forgeの意味がわからなかったのもcondaを敬遠した理由でした。
要はconda公式リポジトリにそのパッケージはないけど、ユーザーグループが管理している conda-forge にはある時にこれを使います。

さて、チートシートに主なコマンドはまとまっているので、詳しいことはそちらを見てもらうとして、自分がいつもやっている手順を紹介します。

最初にインストールしたいパッケージの名前を確認しておきます。


# インストール済みのパッケージの中に、該当のパッケージがないことを確認する。(必要に応じてgrep)
conda list
# 公式リポジトリ内に存在するかどうか確認
conda search [パッケージ名]
# 公式リポジトリにあったらそこからインストール
conda install [パッケージ名]
# 公式リポジトリになかったら conda-forge 内を検索
conda search -c conda-forge [パッケージ名]
# conda-forgeに存在したらそこからインストール
conda install -c conda-forge [パッケージ名]
# バージョンを指定したいときは==の後ろにバージョン番号を指定
conda install [パッケージ名]==[バージョン番号]

また、インストール済みのパッケージのアップデートとアンインストールは次です。


# アップデート
conda update [パッケージ名]
# アンインストール
conda uninstall [パッケージ名]

前置きが長くなって、記事が冗長になってしまったので一旦ここまで。
次はconda-forgeでも見つからないライブラリをPyPIからcondaで入れることのできる、
conda buildや conda skeletonについて紹介したいです。

Python(statsmodels)でインパルス応答関数の計算

思いのほか理論編が長くなってしまいすでにインパルス応答で4記事目ですが、ようやく実装の話です。
グレンジャー因果性と同様に、statsmodelsで学習したVARモデルはインパルス応答関数のメソッドも持っています。

メソッド名は irf です。
statsmodels.tsa.vector_ar.var_model.VARResults.irf

早速やってみましょう。
モデルは、直交化インパルス応答関数の計算例の記事と同じです。
サンプルデータの生成と、パラメーターの推定まで完了しているとします。
今回は次の状態になりました。


# 学習したモデルのパラメーター
result = model.fit(maxlags=1)
print(result.summary())
"""
  Summary of Regression Results   
==================================
Model:                         VAR
Method:                        OLS
Date:           Sun, 23, Feb, 2020
Time:                     21:45:36
--------------------------------------------------------------------
No. of Equations:         2.00000    BIC:                   0.924001
Nobs:                     99.0000    HQIC:                  0.830357
Log likelihood:          -312.903    FPE:                    2.15278
AIC:                     0.766721    Det(Omega_mle):         2.02800
--------------------------------------------------------------------
Results for equation y1
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const         0.177695         0.918241            0.194           0.847
L1.y1         0.553623         0.112916            4.903           0.000
L1.y2         0.169688         0.169945            0.998           0.318
========================================================================

Results for equation y2
========================================================================
           coefficient       std. error           t-stat            prob
------------------------------------------------------------------------
const         1.306800         0.408504            3.199           0.001
L1.y1         0.085109         0.050234            1.694           0.090
L1.y2         0.760450         0.075605           10.058           0.000
========================================================================

Correlation matrix of residuals
            y1        y2
y1    1.000000  0.629788
y2    0.629788  1.000000
"""

サマリーの最後に、撹乱項の相関行列がありますが、分散共分散行列がみたいのでそちらも表示しておきましょう。
sigma_uに入ってます。


# 分散共分散行列
print(result.sigma_u)
"""
[[4.12100608 1.1546159 ]
 [1.1546159  0.81561201]]
"""

さて、まずは非直行化インパルス応答関数を計算して可視化してみましょう。
ちなみに、直行化/非直行化を指定するのは可視化(plot)の引数であり、計算時点で両方とも算出されています。


import matplotlib.pyplot as plt

# 10期先までの非直交化インパルス応答関数のプロット
period = 10
irf = result.irf(period)
irf.plot()
plt.show()

結果がこちらです。

想定通りの結果になっていますね。

続いて、(直行化)インパルス応答関数です。
用語としては、単にインパルス応答関数と言ったら、直行化の方をさすのですが、ライブラリでは明示的に、
orth=Trueを指定しなければなりません。(デフォルトはFalse)

計算は完了しているので、可視化のみです。


irf.plot(orth=True)
plt.show()

結果がこちら。

0期先の値が特徴的で、 y2からy1へは影響していないのに、y1からy2へは影響が発生していますね。
また、もう一点注意しないといけない点があります。
非直行化の方は、y1からy1 および y2からy2への0期先の影響の値が1であることから分かる通り、
1単位のショックを与えた時のインパルス応答関数を計算されていました。
しかし、直行化の方は、そうではなく、1標準偏差のショックを与えた時のインパルス応答関数の値になっています。
(この辺の挙動はドキュメントのどこに書いてあるのだろう?)

結果をプロットするのではなく、値を利用したい場合もあると思います。
その場合、次の2変数にプロットされた値がそれぞれ入っていますのでこれが使えます。


irf.irfs
irf.orth_irfs

直交化インパルス応答関数の計算例

今回は直交化インパルス応答関数を実際に計算してみます。

例として扱うのは、非直交化インパルス応答関数の記事と同じ次の例です。
$$
\left\{\begin{matrix}\
y_{1t}&=& -1+ 0.6 y_{1,t-1}+ 0.3 y_{2,t-1}+\varepsilon_{1t}\\\
y_{2t}&=& 1 + 0.1 y_{1,t-1}+ 0.8 y_{2,t-1}+\varepsilon_{2t}\
\end{matrix}\
\right.\
,\left(\begin{matrix}\varepsilon_{1t}\\\varepsilon_{2t}\end{matrix}\right)\sim W.N.(\boldsymbol{\Sigma})
$$
$$
\boldsymbol{\Sigma} = \left(\begin{matrix}4&1.2\\1.2&1\end{matrix}\right)
$$

さて、まずは$\boldsymbol{\Sigma}$を分解する必要があります。

$$
\begin{align}
\mathbf{A}=
\left(\begin{matrix}1&0\\0.3&1\end{matrix}\right),\\
\mathbf{D}= Var(\mathbf{u}_t)
\left(\begin{matrix}4&0\\0&0.64\end{matrix}\right)
\end{align}
$$
とすると、
$$
\boldsymbol{\Sigma} = \mathbf{ADA^{\top}}
$$
が成立します。
このため、$\boldsymbol{\varepsilon_t}$は次のように、
直交化撹乱項に分解できます。
$$
\left\{\begin{align}\
\varepsilon_{1t}&=u_{1t}\\\
\varepsilon_{2t}&=0.3u_{1t}+u_{2t}\
\end{align}\right.
$$
この撹乱項を用いることで、元のモデルがこのように変わります。

$$
\left\{\begin{matrix}\
y_{1t}&=& -1&+ &0.6 y_{1,t-1}&+ &0.3 y_{2,t-1}&+&u_{1t}\\\
y_{2t}&=& 1 &+ &0.1 y_{1,t-1}&+ &0.8 y_{2,t-1}&+&0.3u_{1t}+u_{2t}\
\end{matrix}\
\right.\
$$
$$
\mathbf{u}_t \sim W.N.(\mathbf{D})
$$
$$
\mathbf{D} = \left(\begin{matrix}4&0\\0&0.64\end{matrix}\right)
$$

この形に変形できれば、非直交化インパルス応答関数と同じ様に計算が出来ます。
例えば、$y_1$に1単位のショックを与えたときの同時点におけるインパルス応答は次になります。
$$
\begin{align}
IRF_{11}(0) = \frac{\partial y_{1t}}{\partial u_{1t}} = 1\\
IRF_{21}(0) = \frac{\partial y_{2t}}{\partial u_{1t}} = 0.3
\end{align}
$$
$IRF_{21}(0)\neq0$となるのが大きな特徴です。

ただし、逆に$y_2$に1単位のショックを与えても、$y_1$には影響がありません。
要するに次のようになります。
$$
\begin{align}
IRF_{12}(0) = \frac{\partial y_{1t}}{\partial u_{2t}} = 0\\
IRF_{22}(0) = \frac{\partial y_{2t}}{\partial u_{2t}} = 1
\end{align}
$$
これが三角分解の仮定により発生している現象であり、変数の並べ方が結果に影響してしまっている点です。
このためインパルス応答関数を使う時は、変数の並べ方に気を使う必要があります。

1期後先以降の値は、非直交化インパルス応答関数のときと同じ漸化式で逐次的に求めることが可能です。

直交化インパルス応答関数

前回の記事: 非直交化インパルス応答関数
の続きです。今回も引用元は沖本本です。

最後に書いたとおり、非直交化インパルス応答関数には、撹乱項の間の相関を無視しているという問題がありました。
これを解決するには、撹乱項を互いに無相関な撹乱項に分解して、それぞれにショックを与えたときの各変数の変化を調べることです。
しかし、一般に撹乱項を互いに無相関な撹乱項に分解することは不可能で、何かしらの仮定が必要になります。
そこで、撹乱項の分散共分散行列$\boldsymbol{\Sigma}$の三角分解を用いて、撹乱項を互いに無相関な撹乱項に分解できると仮定します。

定義 (直交化インパルス応答関数):
撹乱項の分散共分散行列の三角分解を用いて、撹乱項を互いに無相関な撹乱項に分解したとき、
互いに無相関な撹乱項は直交化撹乱項といわれる。また、$y_j$の直交化撹乱項に1単位または1標準偏差のショックを与えた時の、
$y_i$の直交化インパルス応答を時間の関数としてみたものは、$y_j$に対する$y_i$の直交化インパルス応答関数といわれる。

一般的に、インパルス応答(関数)というというと直交化インパルス応答(関数)を指すそうです。
(後日の記事でstatsmodelsによる実装紹介しますが、このモジュールでは非直交化のほうをデフォルトにしてるように見えます。
もしそうだとしたら気をつけて使わないといけないので、記事にする前にもう一回検証します。)

さて、もう少し具体的に定義します。
まず撹乱項 $\boldsymbol{\varepsilon}_t$ の分散共分散行列$\boldsymbol{\Sigma}$は正定値行列であるので、
$\mathbf{A}$を対角成分が1に等しい下三角行列、$\mathbf{D}$を対角行列として、
$$
\boldsymbol{\Sigma} = \mathbf{ADA^{\top}}
$$
と三角分解することが出来ます。
(と、本には書いてありますが、分散共分散行列が正定値なだけでなく、対称行列であることも考慮する必要あると思います。)

このとき、直交化撹乱項 $\mathbf{u}_t$は、
$$
\mathbf{u}_t = \mathbf{A}^{-1}\boldsymbol{\varepsilon}_t
$$
で定義され、各変数固有の変動を表すとみなされます。
なお、
$$
\begin{align}
Var(\mathbf{u}_t) &= Var(\mathbf{A}^{-1}\boldsymbol{\varepsilon}_t) = \mathbf{A}^{-1} Var(\boldsymbol{\varepsilon}_t)(\mathbf{A}^{-1})^{\top} \\
& =\mathbf{A}^{-1}\boldsymbol{\Sigma}(\mathbf{A}^{\top})^{-1}=\mathbf{D}
\end{align}
$$
であることから、$\mathbf{u}_t$が互いに無相関であることが分かります。
また、それぞれの分散も確認できます。

$y_j$のショックに対する$y_i$の(直交化)インパルス応答関数は、$u_{jt}$に1単位のショックを与えたときの、
$y_i$の反応を時間の関数と見たものであり、
$$
IRF_{ij}(k) = \frac{\partial y_{i,t+k}}{\partial u_{jt}}, \ \ \ k=0, 1, 2, \dots.
$$
と表されます。
1単位ではなくて、1標準偏差のショックを与えたときの変数の反応を見たい時は、
$IRF_{ij}(k)$に$\sqrt{Var(u_{jt})}$を掛けることで求まります。(非直交化のときと同じですね)

1標準偏差のショックを与えたときのインパルス応答関数を求める方法として、
三角分解の変わりにコレスキー分解を用いて撹乱項を分解して、
その分解した撹乱項に1単位ショックを与える方法も紹介されていますが、
撹乱項の絶対値が標準偏差倍違うだけの話なので省略します。