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

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

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

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

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

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

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

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

  1. matplotlibのグラフを高解像度で保存する
  2. macにgraphvizをインストールする
  3. pythonで累積和
  4. DataFrameを特定の列の値によって分割する
  5. INSERT文でWITH句を使う
  6. kerasのto_categoricalを使ってみる
  7. numpyのpercentile関数の仕様を確認する
  8. scipyで階層的クラスタリング
  9. matplotlibのデフォルトのフォントを変更する
  10. graphvizで決定木を可視化

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

PrestoのUNNESTを利用した横縦変換

以前Prestoのクエリで縦横変換(pivot)を行う方法を初回しましたが、今回はその逆で横縦変換(unpivot)を紹介します。

参考: PrestoのMap型を使った縦横変換

参考記事と逆の変換やるわけですね。
そのため元のテーブルがこちら。

横持ちのテーブル
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

一番シンプルな書き方は、UNIONを使う方法だと思います。
key の値ごとにvalue を抽出してそれぞの結果を縦に積み上げます。


SELECT
    uid,
    'c1' AS key,
    c1 AS value
FROM
    htable
UNION ALL SELECT
    uid,
    'c2' AS key,
    c2 AS value
FROM
    htable
UNION ALL SELECT
    uid,
    'c3' AS key,
    c3 AS value
FROM
    htable

ただ、この書き方には課題もあって、元の列数が多いとクエリが非常に冗長になります。
そこで、 UNNESTを使った方法を紹介しておきます。

ドキュメントは SELECT の説明のページの途中に UNNESTの章があります。

これを使うと次の様に書けます。


SELECT
    uid,
    t.key,
    t.value
FROM
    htable
CROSS JOIN UNNEST (
  array['c1', 'c2', 'c3'],
  array[c1, c2, c3]
) AS t (key, value)

とてもシンプルに書けました。

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(黒い点)が配置できています。
ただ、全部相異なることを確認するのは少し手間ですね。

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

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で表したスートの2文字で表したもの)を強い順にソートするものです。


# カードの強さを返す関数
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を使って、アルファベットを小文字に統一してソートするなど、
使い方は色々ありそうです。

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で十分なことが多いと思います。)

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

先ほどの 6種類のカーネル + 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が頑張ってますね。
全データを学習に使えるのも大きいかもしれません。

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のコードもかなり参考になりました。

さて、推定したバンド幅はその2乗の値が covariance という属性に格納されています。(1次元データの場合。)
これは、 学習したデータの不変分散に、 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乗を掛けてみると、先ほどのバンド幅の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乗を掛けてみると、先ほどのバンド幅の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を使った実装を紹介する予定です。

SciPyによるカーネル密度推定

前回の記事に続いてカーネル密度推定の話です。
参考: カーネル密度推定法の紹介

前の記事では数式に沿って実装するため、numpyで書きましたが、普段の利用では、SciPyにクラスが実装されているのでこれを使えます。

参考: scipy.stats.gaussian_kde
クラス名から分かる通り、カーネルはガウスカーネルがあらかじめ指定されており、他の関数は使えません。
実際のところ、ガウスカーネル以外が使えなくて困ることは滅多にないでしょう。

使い方ですが、 インスタンスを作るときに、データセットを渡すだけで完了します。
とてもありがたいことに、バンド幅の指定が要りません。
デフォルトでは Scottのルールと呼ばれる方法で決定されます。
bw_methodという引数で指定することで、もっと直接的な指定や、Silvermanのルールを使うことも可能です。
ScottのルールとSilvermanのルールがそれぞれ具体的にどの様にバンド幅を決めているのかは後日のブログで紹介するとして、
とりあえず、正規分布から生成したデータで使ってみましょう。


from scipy.stats import gaussian_kde
from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt

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

# グラフ描写のため真の確率密度関数の値を取得しておく
xticks = np.linspace(-5, 9, 141)
pdf_values = norm(loc=2, scale=2).pdf(xticks)

# カーネル密度推定
kde = gaussian_kde(data)
# evaluate には pdfというエイリアスがあるのでそちらでも可
estimated_value = kde.evaluate(xticks)

# 可視化
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)
ax.plot(xticks, pdf_values, label="真の分布")
ax.plot(xticks, estimated_value, label="推定値",  linestyle="--")
# 標本データ
for x in data:
    ax.vlines(x, ymin=0, ymax=0.01, alpha=0.5)
ax.set_ylim((0, 0.25))
ax.legend()
plt.show()

出力結果がこちらです。

とても手軽に良い推定が得られていますね。

ちなみに、推定してくれたバンド幅ですが、covariance という属性にその情報を持っています。


# バンド幅の確認 (hの2乗なので注意)
print(kde.covariance)

# [[0.79905207]]

多次元にも対応するため、分散共分散行列の形でデータを持てる様になっています。
今回は1次元データなので、ちょっと面倒ですが、無駄に2次元配列になってる中から要素を持ってくる必要があります。
varianceなので、この値は分散、前回の記事のバンド幅は標準偏差なので、平方根を取る必要があることも注意です。

せっかく推定してくれたバンド幅が取れるので、前回スクラッチで描いたコードの推定値と結果が一致することも見ておきましょう。
前回の記事のコピーだと、関数名が重複してしまって面倒なので適当に、my_kdeにリネームしました。


# ガウスカーネル
def gaussian_kernel(u):
    return(np.exp(-(u**2)/2)/np.sqrt(2*np.pi))


# ガウスカーネルを用いたカーネル密度推定
def my_kde(x, h, dataset):
    return sum([gaussian_kernel((x-d)/h) for d in dataset])/(len(dataset)*h)


# バンド幅の取得
h = np.sqrt(kde.covariance[0][0])
# 可視化
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)
ax.plot(xticks, pdf_values, label="真の分布")
ax.plot(xticks, kde.evaluate(xticks), label="scipy", linestyle="--")
ax.plot(xticks, my_kde(xticks, h, data), label="スクラッチ", linestyle=":")
ax.legend()
plt.show()

SciPy で出した結果と、スクラッチで書いてるやつの結果が(一致してることの確認なので)重なってしまっていて、見にくいですが、
一応線の種類を変えてあるので、よく見るとしっかり同じところに線が引かれているのが分かります。

カーネル密度推定法の紹介

今回からしばらくカーネル密度推定法について書いていきます。

参考文献ですが「パターン認識と機械学習 上 (PRML 上)」の 119ページ、 2.5.1 カーネル密度推定法 あたりがわかりやすいのではないでしょうか。
ただし、この本だと最後にパラメーター$h$の決め方が重要だ、ってところまで書いて、$h$の決め方を紹介せずに次の話題に移っているので、
これからの数回の記事の中でそこまで説明できたらいいなと思っています。

さて、カーネル密度推定法を大雑把に言うと、
「確率密度が全然わからない確率分布から標本が得られたときに、その標本からできるだけいい感じに元の確率密度を推定する方法」です。
元の確率密度が正規分布だとかポアソン分布だとかだけでもわかっていたら、その期待値や分散などのパラメーターだけ推定すればいい(これをパラメトリックという)のですが、
それさえもわからないときに使われる手法であり、ノンパラメトリックなアプローチと呼ばれます。

できるだけ一般的な形でカーネル密度推定法の方法を書くと次のような形になります。

ある$D$次元のユークリッド空間の未知の確率密度$p(\mathbf{x})$から、観測値の集合$\{\mathbf{x}_n\}\ \ (n=1,2, \cdots, N)$が得られたとします。

そして、次の性質を満たすカーネル関数$k(\mathbf{u})$を選びます。
$$
\begin{align}
k(\mathbf{u}) &\geq 0,\\
\int k(\mathbf{u}) d\mathbf{u} &= 1
\end{align}
$$
このとき、$p(\mathbf{x})$を次で推定するのがカーネル密度推定法です。($h$はパラメーター。)
$$
p(\mathbf{x}) = \frac{1}{N}\sum_{n=1}^{N}\frac{1}{h^D}k(\frac{\mathbf{x}-\mathbf{x}_n}{h}).
$$

カーネル密度関数は要は確率密度間数ですね。理論上は何でもいいそうですが、現実的には原点周りで値が大きく、そこから外れると0に近く(もしくは等しく)なる関数が好都合です。
正規分布の確率密度関数(ガウスカーネル)が頻繁に用いられます。

上記の式が導出されるまでの話はPRMLに書いてあります。
その部分を省略してしまったので代わりに上の定義通り計算して、未知の確率分布が近似できている様子を確認しておきましょう。
できるだけライブラリ使わないでやります。(ただ、正解データの分布作成でscipyは使っていますし、numpyやmatplotlibは例外とします。)
カーネル関数にはガウスカーネルと、単位超立方体内で1、その外で0を取る関数(パルツェン窓というらしい)の2種類でやってみます。

元となる分布は正規分布2個を足した、2つ山のある分布で実験します。($D$は1です)

最初に正解になる分布関数の定義と、そこからのデータサンプリングを済ませておきます。


from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt


# 推定したい分布の真の確率密度関数
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)

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

そして、それぞれのカーネル関数と、それを使ったカーネル密度関数を定義しておきます。


# パルツェン窓
def parzen_window(u):
    if np.abs(u) <= 0.5:
        return 1
    else:
        return 0


# ガウスカーネル
def gaussian_kernel(u):
    return(np.exp(-(u**2)/2)/np.sqrt(2*np.pi))


# パルツェン窓を用いたカーネル密度推定
def parzen_kde(x, h, dataset):
    return sum([parzen_window((x-d)/h) for d in dataset])/(len(dataset)*h)


# ガウスカーネルを用いたカーネル密度推定
def gaussian_kde(x, h, dataset):
    return sum([gaussian_kernel((x-d)/h) for d in dataset])/(len(dataset)*h)

あとは可視化です。 今回は パラメーター$h$は$0.5, 1, 2$の3つ試しました。


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

# 3種類のhで実験
for i, h in enumerate([0.5, 1, 2]):
    parzen_kde_values = [parzen_kde(x, h, X) for x in xticks]
    gaussian_kde_values = [gaussian_kde(x, h, X) for x in xticks]

    ax = fig.add_subplot(3, 2, i*2+1, facecolor="w")
    ax.set_title(f"パルツェン窓, h={h}")
    ax.plot(xticks, pdf_values, label="真の分布")
    ax.plot(xticks, parzen_kde_values, label="パルツェン窓")
    # 標本データ
    for x in X:
        ax.vlines(x, ymin=0, ymax=0.01)
    ax.set_ylim((0, 0.37))
    ax.legend()

    ax = fig.add_subplot(3, 2, i*2+2, facecolor="w")
    ax.set_title(f"ガウスカーネル, h={h}")
    ax.plot(xticks, pdf_values, label="真の分布")
    ax.plot(xticks, gaussian_kde_values, label="ガウスカーネル")
    # 標本データ
    for x in X:
        ax.vlines(x, ymin=0, ymax=0.01)
    ax.set_ylim((0, 0.37))
    ax.legend()

plt.show()

出力結果がこれです。

青い線が真の分布、黒くて短い縦線たちが推定に使った標本データ、
オレンジの線が推定した分布です。

「元の確率密度関数の形が全くわからない」という前提から出発している割にはうまく推定できているのもあると言えるのではないでしょうか。
また、カーネル密度関数の選び方、$h$の決め方がそれぞれ重要であることも分かりますね。

やってることといえば、得られた標本の付近で大きめの値を取る関数を全ての標本について計算してその平均を取っているだけなのですが、シンプルな割に強力です。