カーネル密度関数を書く

1 やりたいこと

Jupyterでカーネル密度関数を生成して、ビジュアライズ化する。
カーネル密度関数の解説はこちら

2 コード

gaussian_kde関数を使うのと、seabonのkdeplot関数を使うのとでカーネル密度関数は若干形状が異なるようだ。

/home/sampletest/sample.py

from numpy.random import randn
import seaborn as sns
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde
sns.set_style("whitegrid")
from scipy.integrate import cumtrapz #pdfを全区間で積分するためのライブラリ



# datasetを与え、カーネル密度関数を発生させる。また、カーネル密度関数の累積度数を計算し、
# 累積値down,upとなるX値を返す。

def pdf_kernel(dataset,down,up):

    d_max=np.max(dataset)+0.2
    d_min=np.min(dataset)-0.2

    d_kernel = gaussian_kde(dataset) #カーネル密度推定関数を生成する。バンド幅は自動計算されている。

    #積分を行う範囲を指定する：d_minからd_maxの範囲で積分を行う。
    d_xs = np.linspace(d_min, d_max, num=1000)

    #カーネル密度関数の入力(d_xs)と出力(d_ys)を定義している。
    d_ys = d_kernel(d_xs)

    #累積分布関数をd_xsの範囲で積分する。
    d_integral = cumtrapz(d_ys, d_xs)

    #cdf(x) = 0.03となるxを求める。d_integral配列の中で0.03に最も近くなる数値が配列中の何番目に当たるかを算出している。
    idx_d= np.searchsorted(d_integral, down)

    #cdf(x) = 0.9となるxを求める。d_integral配列の中で0.9に最も近くなる数値が配列中の何番目に当たるかを算出している。
    idx_u = np.searchsorted(d_integral, up)


    #グラフ表示を行っている。
    ax=plt.plot(d_xs, d_ys, label="KDE")
    plt.xlim(d_min-1, d_max+1) 

    #累積値5%の範囲を表示
    plt.fill_between(d_xs[:idx_d], 0, d_ys[:idx_d], facecolor="r", alpha=0.5)


    #累積値90%の範囲を表示
    plt.fill_between(d_xs[idx_u:], 0, d_ys[idx_u:], facecolor="r", alpha=0.5)

    #凡例を右上に表示
    plt.legend(loc="upper right")

    pdf_val={"down":d_xs[idx_d],"up":d_xs[idx_u]}
    return pdf_val



# メインプログラムは以下から開始
dataset = randn(50) #一様分布に従う乱数を50個発生させる
down=0.05 
up=0.9


# datasetを与え、カーネル密度関数を発生させる。また、カーネル密度関数の累積度数を計算し、
# 累積値down,upとなるX値を返す。
val=pdf_kernel(dataset,down,up) 

# seabornを使い、カーネル密度関数を書く。
sns.kdeplot(dataset,label="from seaborn")

実行結果
#累積値0～down(図示例では5%)の範囲と、累積値up(図示例では90%)～100%の範囲を表示している。