分析ノート

2020-07-132020-07-14

jupyter notebookの拡張機能を導入する

職場のMacでは以前から導入していたのですが、私物のMacにも入れたのでメモしておきます。

そのままでも十分使いやすい jupyter notebookですが、
jupyter_contrib_nbextensions　というライブラリを入れるといろいろな拡張機能が使える様になり、一層便利になります。
自動で pep8 を満たす様に整形してくれたり、入力セルを隠せたり、補完機能が強化されたり、
テーブルがきれいに表示されるなどの設定ができます。

インストールはこちらのページにもある通り、
pipでも condaでも入れることができます。自分はnotebook自体がcondaで導入したものなので、
condaで入れています。
これを jupyter notebookが起動してない時に実行し、その後、notebookを起動します。


# pipで入れる場合のコマンド
# $ pip install jupyter_contrib_nbextensions
# condaで入れる場合のコマンド
$ conda install -c conda-forge jupyter_contrib_nbextensions

すると、 jupyterの画面に、 Nbextensions というタブが生成されます。

チェックボックスのチェックを外すと設定を変えられるようになり、
あとは好みの設定を選んで入れていくだけです。

実際に notebookを触りながら色々試すと楽しいと思います。

(補足) サイトによっては、利用するため(タブを表示する)にはアクティベーションの操作が必要といったことが書いてありますが、
僕の環境ではインストールしたらタブが勝手に出てきました。
バージョンによって違うのかもしれません。

2020-07-092020-07-04

Googleアナリティクスでアクセスの増加量が大きいページを探す

2020年の5月4日ごろから、Googleの検索アルゴリズムのコアップデートが行われ、このブログの訪問者数にも大きな影響が出ています。
アップデートの前後3週間の、自然検索による流入数をグラフにしたのが下図です。
明らかに流入が増えていますね。

さて、もっと詳しく、具体的にどのページの流入が増えて、どのページの流入は減ったのか調べるため、方法を調べたので紹介します。
デフォルトだと、多い順になっており、逆順になれべても絶対値の少ない順になってしまって、変化量順にはならないですよね。

変化量でソートしたい場合は、「並べ替えの種類:」という項目を「変化量」に設定する必要があります。
これで、アクセスの増加が激しいページを特定できました。
この期間はゴールデンウィークなどの要因もあるので、変化量は注意しいてみないといけないのですが、
なんとなくテクニカルな記事が恩恵を受けてるように感じます。

2020-07-022021-06-27

2020年上半期(1月~6月)によく読まれた記事

早いもので2020年が半分終わってしまいました。
今年はブログ更新頻度を落としているのもあり、四半期でのランキング発表をやめているので、
この辺で半年間に読まれた記事のランキングを出したいと思います。

参考ですが、昨年1年間のランキングはこちらです。
参考: 2019年のまとめ

では早速ランキング発表です。
集計期間は 2020年1月から6月まで。pvで並べています。

データサイエンスというより、プログラミングのちょっとしたTips的な記事の方がよくpvを集めていますね。
今年はネットワーク解析/グラフ理論の記事も頑張って書いたので今後はそれらのランクインも期待したいです。

2020-06-292020-06-28

横持ちのテーブル
uid	c1	c2	c3
101	11	12	13
102	21	22	23

縦持ちのテーブル (vtable)
uid	key	value
101	c1	11
101	c2	12
101	c3	13
102	c1	21
102	c2	22
102	c3	23

NetworkXを使ってエイト・クイーンパズルを解く

“Python言語によるビジネスアナリティクス”と言う本の中に、NetworkXについての記述があるのですが、
その演習問題(問題43)として8-クイーン問題が出題されています。
面白そうだったので挑戦してみたところ、なんとか解けました。
この問題がネットワーク解析に関係あるとは意外ですね。

まず、8-クイーン問題については、Wikipediaの説明が分かりやすいと思います。
エイト・クイーン – Wikipedia

簡単に言うと、8×8の盤面状にチェスのクイーン(将棋の飛車と角を足した動きができる最強のコマ)をお互いに取り合えない位置に並べるパズルです。
正解は92通り(盤面の反転や回転を除いた本質的なものは12通り)あるそうです。

さて、これを解くために、次のアプローチを考えました。

1. 8×8の各マスをノードとするグラフを構築する。
2. クイーンがお互いに取り合えるマス同士をエッジで結ぶ。
3. そのグラフの補グラフを取得することで、クイーンがお互いに取り合えないマス通しが結ばれたグラフを作る。
4. ノード数が8のクリーク(その中に含まれる全てのノードが結ばれている部分グラフ)を探索する。

最初からお互いに取り合えないマスの間を結ぶのではなく、2.と3.に分けてるのはその方が実装がシンプルになるからです。


import networkx as nx

# 盤面の一辺のマス数
bord_size = 8

# グラフ生成
G = nx.Graph()

# ノードを作成する
for i in range(bord_size):
    for j in range(bord_size):
        G.add_node((i, j))

# 同じ行か同じ列に並ぶノードをエッジで結ぶ
for i in range(bord_size):
    for m in range(bord_size):
        for n in range(m+1, bord_size):
            G.add_edge((i, m), (i, n))
            G.add_edge((m, i), (n, i))

# 同じ斜め線状に並ぶノードをエッジで結ぶ
for n1 in G.nodes:
    for n2 in G.nodes:
        if n1 == n2:
            continue
        elif n1[0]+n1[1] == n2[0]+n2[1]:
            G.add_edge(n1, n2)
        elif n1[0]-n1[1] == n2[0]-n2[1]:
            G.add_edge(n1, n2)


# 補グラフを得る (お互いに取り合わない辺が結ばれている)
G_complement = nx.complement(G)

# サイズが変の数に等しいクリークの一覧を得る
answers = [
        clieque for clieque in nx.find_cliques(G_complement)
        if len(clieque) == bord_size
    ]

# 得られた解の個数
print(len(answers))
# 92

こうして、92個の解が得られました。
あとは実際にこれを可視化してみましょう。


import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(12, 30), facecolor="w")

for i, answer in enumerate(answers, start=1):
    bord = np.zeros([bord_size, bord_size])
    ax = fig.add_subplot(16, 6, i)

    for cell in answer:
        bord[cell] = 1

    ax.imshow(bord, cmap="gray_r")

plt.show()

出力がこちらです。　一個一個みていくと確かに縦横斜めにダブらない様にQueen(黒い点)が配置できています。
ただ、全部相異なることを確認するのは少し手間ですね。

もともと出題されていた本には解答が載っていないので、
もしかしたらこの記事のコードでは相当非効率なことをやっているかもしれませんが、
おそらくこんなイメージで解くのだと思います。

2020-06-222020-06-21

Pythonのlistをsortするときに指定するkey引数について

先日とある目的で、配列を要素に持つ配列を、各配列の3番目の値を使ってソートするコードが必要になりました。
最初はそのキーになる要素だけ取り出して、 argsortしてどうにか並べ替えようとしたのですが、実はkeyって引数を使えば綺麗に実装できることがわかりました。

ドキュメントは組み込み型のlistのsortメソッドの部分が該当します。

key って引数に単一の引数を取る関数を渡せば、listの各要素にその関数が適用され、戻り値によってソートされる様です。
なので、配列の配列を3番目の要素でソートするには、lambda 式あたりでそう言う関数を作ってあげれば実現できます。
(インデックスは0始まりなので、実際にはインデックス2を取得します。)


# ソートされるデータ作成
list_0 = [
        [9, 8, 5],
        [0, 8, 3],
        [1, 6, 5],
        [9, 0, 0],
        [4, 9, 3],
        [1, 4, 8],
        [4, 0, 6],
        [0, 3, 5],
        [1, 3, 1],
        [5, 2, 7],
    ]

# 3番目(index 2)の要素でソート
list_0.sort(key=lambda x: x[2])
print(list_0)
# [[9, 0, 0], [1, 3, 1], [0, 8, 3], [4, 9, 3], [9, 8, 5], [1, 6, 5], [0, 3, 5], [4, 0, 6], [5, 2, 7], [1, 4, 8]]

これを応用すれば、辞書の配列を特定のkeyに対応する値でソートするといったことも簡単にできます。
(lambda式に辞書オブジェクトのgetメソッド渡すだけですね。)

もちろん、lambda 式以外にも事前に定義した関数を渡すこともできます。
次の例は、トランプのカード (2〜9とTJQKAで表したランクと、cdhsで表したスートの２文字で表したもの)を強い順にソートするものです。


# カードの強さを返す関数
def card_rank(card):
    rank_dict = {
            "T": 10,
            "J": 11,
            "Q": 12,
            "K": 13,
            "A": 14,
        }
    suit_dict = {
            "s": 3,
            "h": 2,
            "d": 1,
            "c": 0,
        }
    return int(rank_dict.get(card[0], card[0]))*4 + suit_dict[card[1]]


card_list = ["Ks", "7d", "Ah",  "3c", "2s", "Td", "Kc"]

# ソート前の並び
print(card_list)
# ['Ks', '7d', 'Ah', '3c', '2s', 'Td', 'Kc']

# 強い順にソート
card_list.sort(key=card_rank, reverse=True)
print(card_list)
# ['Ah', 'Ks', 'Kc', 'Td', '7d', '3c', '2s']

reverse = True を指定しているのはソートを降順にするためです。

このほかにも str.lowerを使って、アルファベットを小文字に統一してソートするなど、
使い方は色々ありそうです。

2020-06-182020-06-17

statsmodels でカーネル密度推定

今回は statsmodelsを使ったカーネル密度推定を紹介します。

ドキュメントはこちらです: statsmodels.nonparametric.kde.KDEUnivariate

statsmodelsのいつもの作法で、インスタンスを作る時にデータを渡し、fitメソッドを呼び出せば完了です。
(SciPyではfitは必要なかったので注意)

fit する時に、カーネル関数の種類やバンド幅の決定方法を指定できます。
詳しくは fit メソッドのドキュメントを確認してください。

カーネル関数が 7種類指定できたり、バンド幅の指定方法が少し違ったりします。

ではやっていきましょう。比較のため、 SciPyで推定した結果も載せていきます。
(SciPyのバンド幅指定のルールがおかしいと言う話もあるので。参考記事はこちら。)
まずはデータ生成です。


import numpy as np
import statsmodels.api  as sm
import matplotlib.pyplot as plt
from scipy.stats import norm, gaussian_kde

# 推定したい分布の真の確率密度関数
def p_pdf(x):
    return (norm.pdf(x, loc=-4, scale=2)+norm.pdf(x, loc=2, scale=1))/2


# その分布からのサンプリング
def p_rvs(size=1):
    result = []
    for _ in range(size):
        if norm.rvs() < 0:
            result.append(norm.rvs(loc=-4, scale=2))
        else:
            result.append(norm.rvs(loc=2, scale=1))
    return np.array(result)


# グラフ描写のため真の確率密度関数の値を取得しておく
xticks = np.linspace(-10, 5, 151)
pdf_values = p_pdf(xticks)

# 100件のデータをサンプリング
data = p_rvs(100)

さて、ここから SciPyと statsmodelsでそれぞれのScottのルールとSilvermanのルールでカーネル密度推定して結果を可視化します。


# SciPyのカーネル密度推定
kde_scipy_scott = gaussian_kde(data, bw_method="scott")
kde_scipy_silverman = gaussian_kde(data, bw_method="silverman")

# statusmodels のカーネル密度推定
kde_sm_scott = sm.nonparametric.KDEUnivariate(data)
kde_sm_scott.fit(kernel="gau", bw="scott")
kde_sm_silverman = sm.nonparametric.KDEUnivariate(data)
kde_sm_silverman.fit(kernel="gau", bw="silverman")

fig = plt.figure(figsize=(7, 7), facecolor="w")
ax = fig.add_subplot(1, 1, 1, title="カーネル密度推定の結果")
ax.plot(xticks, pdf_values, label="真の分布")
ax.plot(xticks, kde_scipy_scott.evaluate(xticks), label="SciPy - Scott")
ax.plot(xticks, kde_scipy_silverman.evaluate(xticks), alpha=0.5, label="SciPy - Silverman")
ax.plot(xticks, kde_sm_scott.evaluate(xticks), linestyle=":", label="statsmodels - Scott")
ax.plot(xticks, kde_sm_silverman.evaluate(xticks), label="statsmodels - Silverman")
ax.hist(data, bins=20, alpha=0.1, density=True, label="標本")
ax.legend()
plt.show()

出力がこちらです。

statsmodels で Silverman のルールを使ったのが一番良さそうですね。
ただ、これは常にそうなるわkではなく、元の分布の形や標本サイズによって何が最適かは変わってきます。

SciPyのSilvermanのルールと、statsmodelsのScottのルールが重なってしまうのは、先日の記事で書いた通り、SciPyの不備と思われます。

さて、こうしてみると、 statsmodels が一番良さそうに見えますが、当然これにも欠点があります。
SciPyやscikit-learnのカーネル密度推定では多次元のデータも扱えるのですが、
statsmodelsでは1次元のデータしか扱えません。

カーネル密度推定に関しては SciPy, scikit-learn statusmodels　それぞれにメリットデメリットあるので、
目的に応じて使い分けていくのが良さそうです。
(これにそこまで高い精度が求められることも少ないので一番手軽に書ける SciPyで十分なことが多いと思います。)

2020-06-152020-06-13

scikit-learnでカーネル密度推定

今回もカーネル密度推定です。前回の記事の予告通り、今回はscikit-learnを使います。

ドキュメントはこちらです: sklearn.neighbors.KernelDensity

scikit-learnの他のモデルと同様に、import してインスタンス作って、fitしたら学習完了、となるのですが、その後の推定に注意することがあります。
version 0.23.1 時点では、他のモデルにある、predict やtransform がありません。
これは代わりに、score_samples を使います。そして戻り値は、
Evaluate the log density model on the data.
とある様に、その点での密度の”対数”です。なので、SciPyの時みたいに密度関数を得たかったら指数関数で変換する必要があります。

さて、scikit-learn でカーネル密度推定をやるメリットですが、ガウスカーネルだけでなく、全部で6種類のカーネル関数が使えます。

1点だけのデータにたいしてカーネル密度推定を実施すると、そのカーネルの形がそのまま残るのでそれを利用してカーネル関数を可視化しました。


import numpy as np
import matplotlib.pyplot as plt
from sklearn.neighbors import KernelDensity

kernel_names = [
        "gaussian",
        "tophat",
        "epanechnikov",
        "exponential",
        "linear",
        "cosine",
    ]

# x軸のメモリようの配列
xticks = np.linspace(-3, 3, 61)

fig = plt.figure(figsize=(10, 7), facecolor="w")
for i, kernel in enumerate(kernel_names, start=1):
    # 指定したカーネルでモデル作成
    kde_model = KernelDensity(bandwidth=1, kernel=kernel)
    # 1点だけのデータで学習
    kde_model.fit([[0]])
    ax = fig.add_subplot(2, 3, i, title=kernel)
    ax.set_ylim([0, 1.1])
    ax.plot(xticks, np.exp(kde_model.score_samples(xticks.reshape(-1, 1))))

plt.show()

結果がこちら。それぞれ個性ありますね。

さて、カーネル密度推定はscikit-learnで行う場合もバンド幅の決定の問題があります。
そして、SciPyの様に自動的には推定してくれません。

そこで、scikit-learnでは、他のモデルと同じ様に、最適なパラメーターを探索することになります。
Score というメソッドで、与えられたデータに対する対数尤度の合計値が算出できるので、これを使ってグリッドサーチすると良さそうです。
他のどんなモデルもそうなのですが、訓練データと評価データはかならず分ける必要があります。
特にカーネル密度推定では、バンド幅が狭ければ狭いほど訓練データに対するスコアは高くなり、極端に過学習します。
今回は久々にグリッドサーチでもやりましょうか。
また、正解の分布はベータ分布を使います。 (正規分布だとガウスカーネルが有利すぎる気がしたので)

先ほどの６種類のカーネル + SciPyでやってみました。


from scipy.stats import beta
from sklearn.model_selection import GridSearchCV

# グリッドサーチでバンド幅決めてベストのモデルを返す関数
def search_bandwidth(data, kernel="gaussian", cv=5):
    bandwidth = np.linspace(0.1, 2, 190)
    gs_model = GridSearchCV(
                KernelDensity(kernel=kernel),
                {'bandwidth': bandwidth},
                cv=cv
            )
    gs_model.fit(data.reshape(-1, 1))
    return gs_model.best_estimator_

# 正解の分布
frozen = beta(a=3, b=2, loc=-1, scale=2)
# データ生成
data = frozen.rvs(200)

# x軸のメモリようの配列
xticks = np.linspace(-1.1, 1.1, 221)


fig = plt.figure(figsize=(10, 11), facecolor="w")
for i, kernel in enumerate(kernel_names, start=1):
    # 指定したカーネルでモデル作成
    kde_model = search_bandwidth(data, kernel=kernel)
    bw_str = str(kde_model.bandwidth.round(2))
    ax = fig.add_subplot(3, 3, i, title=kernel + " h=" + bw_str)
    # ax.set_ylim([0, 1.1])
    ax.plot(xticks, np.exp(kde_model.score_samples(xticks.reshape(-1, 1))))
    ax.plot(xticks, frozen.pdf(xticks), linestyle="--")

# おまけ。scipy
kde = gaussian_kde(data)
ax = fig.add_subplot(3, 3, 7, title="SciPi")
ax.plot(xticks,  kde.evaluate(xticks))
ax.plot(xticks, frozen.pdf(xticks), linestyle="--")

plt.show()

結果がこちら。一番手軽に使える SciPyが頑張ってますね。
全データを学習に使えるのも大きいかもしれません。

2020-06-112020-06-13

SciPyのガウシアンカーネル密度推定に実装されているバンド幅の算出方法について

僕の環境の SciPyは version 1.4.1 なのでこの記事ではそれを想定して読んでください。(将来のバージョンでは修正されるかも。)

さて、カーネル密度推定においてバンド幅(パラメーターh)を決めるのに使われる、
Scottのルールと、Silvermanのルールについて紹介しようと思っていたのですが、
SciPyの実装を精査していく中でちょっと挙動がおかしいことに気付きました。

それぞれのルールの原典がどっちも無料では見れないのでそれを確認できないのですが、
SciPyで実装されているSilvermanのルールが、どうやら他で言うScottのルールのようです。
SciPyで実装されているScottのルールは、これ一体何なんでしょう。
他で言われているSilvermanのルールはもっと狭いバンド幅で推定するはずです。

Scottのルールの原典:
D.W. Scott, “Multivariate Density Estimation: Theory, Practice, and Visualization”, John Wiley & Sons, New York, Chicester, 1992.
Silvermanのルールの原典:
B.W. Silverman, “Density Estimation for Statistics and Data Analysis”, Vol. 26, Monographs on Statistics and Applied Probability, Chapman and Hall, London, 1986.

本当は正しい(と言うか、この二つのルールも経験則的なルールであって、それを正しいと呼ぶかは疑問符が付きますが)それぞれの発案者が提案したルールを紹介したいのですが、
普段の利用を考えると、scipy.stats.gaussian_kde で、bw_method　に指定した文字列によって挙動がどう変わるのかを正確に理解している方が役に立つと思います。
そのため、この記事のタイトル通り、発案者が何を提唱したかは一旦脇に置いといて、文字列として、”scott”と”silverman”を指定したらプログラムがどう計算するかを紹介します。

ドキュメントはこちらです: scipy.stats.gaussian_kde
GitHubのコードもかなり参考になりました。

さて、推定したバンド幅はその２乗の値が covariance という属性に格納されています。(１次元データの場合。)
これは、学習したデータの不変分散に、 factor という属性の”2乗を”掛けることで計算される実装になっています。
バンド幅hで考えると、不変分散の平方根に factor を掛けたものになりますね。

で、その factor の計算方法が、指定した “scott”と”silverman”で変わります。
データの件数を$n$、データの次元を$d$とすると、
“scott”を指定したときは、
$$
factor = n^{\frac{-1}{d+4}}
$$
です。1次元の場合は$d=1$なので$n^{-0.2}$ですね。

“silverman”を指定したときは、
$$
factor = \left(\frac{d+2}{4}*n\right)^{\frac{-1}{d+4}}
$$
になります。
$d=1$の時は、$(3/4)^{-0.2}=1.05922\cdots$なので、だいたい、$1.06n^{-0.2}$くらいになります。
他のサイトの情報を見ていると、 0.9とか1.06とかの定数が出てくるのですが、そのうち1.06がどこから登場するのかはこれで分かりますね。
ちなみに statsmodelsのカーネル密度推定では 1.059がハードコーディングされています。

で、この 1.06は本当は Scott のルールで使われる定数のはずなのですが、なぜかSciPyでは”silverman”を指定したときに登場します。
1.06の代わりに0.9くらいの値を使うのがSilvermanのルールのはずなのですが。

ただ、以上で、SciPyがどの様にバンド幅を決定しているかはわかったので、その通りの挙動をしていることをコードで確認しておきましょう。
とりあえずデータを準備します。


from scipy.stats import gaussian_kde
from scipy.stats import norm
import numpy as np

# 50件のデータをサンプリング
data = norm(loc=2, scale=2).rvs(50)

print("標本の不偏分散: ", np.var(data, ddof=1))
# 標本の不偏分散:  4.2218467749903645

まずは “scott” の確認です。(初期値なので、コード中に”scott”の文字列は登場しません。)


# Scottのルールでカーネル密度推定
kde = gaussian_kde(data)

# バンド幅 hの2乗を確認
print("バンド幅 h^2: ", kde.covariance[0, 0])
# バンド幅 h^2:  0.8829059945819668

# Scottの因子が、データ件数の -0.2乗に等しいことを確認
print(kde.scotts_factor())
# 0.45730505192732634
print(len(data) ** (-0.2))
# 0.45730505192732634

# データの普遍分散に Scottの因子の2乗を掛けてみると、先ほどのバンド幅の２乗に等しい
print((kde.scotts_factor()**2) * np.var(data, ddof=1))
# 0.8829059945819668

以上の通り、ドキュメント通りの実装になっていました。

次が bw_method=”silverman” の場合です。ごくわずか謎の端数が出てますが、こちらもドキュメント通りの実装であることが確認できます。


# カーネル密度推定 (silvermanのルール)
kde = gaussian_kde(data, bw_method="silverman")

# バンド幅 hの2乗を確認
print("バンド幅 h^2: ", kde.covariance[0, 0])
# バンド幅 h^2:  0.9905809235665322

# Silvermanの因子が、定義通りの計算結果に等しいことを確認
print(kde.silverman_factor())
# 0.48438841363348917
print((kde.n * (kde.d + 2) / 4.)**(-1. / (kde.d + 4)))
# 0.4843884136334891

# データの普遍分散に Silvermanの因子の2乗を掛けてみると、先ほどのバンド幅の２乗に等しい
print((kde.silverman_factor()**2) * np.var(data, ddof=1))
# 0.9905809235665322

ちなみに、 statsmodelsのドキュメントを見てみましょう。
statsmodels.nonparametric.kde.KDEUnivariate.fit

The bandwidth to use. Choices are:
“scott” – 1.059 * A * nobs ** (-1/5.), where A is min(std(X),IQR/1.34)
“silverman” – .9 * A * nobs ** (-1/5.), where A is min(std(X),IQR/1.34)
“normal_reference” – C * A * nobs ** (-1/5.), where C is calculated from the kernel. Equivalent (up to 2 dp) to the “scott” bandwidth for gaussian kernels. See bandwidths.py
If a float is given, it is the bandwidth.

この様に、 “scott” のほうが、 $1.059*n^{(-0.2)}$ を掛けています。SciPyでは”silverman”の時に登場するものですね。
また、データの標準偏差よりも、IQR (これは第3四分位数から第1四分位数を引いた値です)を1.34で割ったものが小さかったらそっちを使うと言うルールも実装されています。
SciPyのコードを見るとこれも入ってないです。

ちょっと色々疑わしいところはあるのですが、ぶっちゃけた話、ヒストグラムなどを使って可視化する時におまけにつけることがあるくらいのものなので、
あんまり気にしないで使って行こうかなと思っています。

ただ、statsmodelsの方が正しい値出してくれそうなのでこれもそのうちコードを載せて紹介しましょう。
次の記事ではガウスカーネル以外を使ったカーネル密度推定を先にやりたいので、Scikit-learnを使った実装を紹介する予定です。

投稿