「カーネル密度推定」について学ぶ(統計学 / 教師なし学習 / 密度推定)

2023年5月13日

この記事では、統計学を初めて学ぶ筆者が、「教師なし学習」における「密度推定」の手法である「カーネル密度推定」について学んだ内容について記載しています。

学習には、Wikipediaの「カーネル密度推定」の記事を参考にし、Pythonのプログラミングにも触れ、理解を深めました。

プログラミングには、機械学習ライブラリのscikit-learnを使用しました。

この記事は、他の人が参考にできるよう、わかりやすく書くことを心がけました。

教師なし学習

教師なし学習は、データからパターンや構造を見つけ出すための機械学習の手法です。データにはラベルや目的変数がなく、機械自身が学習し、データの背後にある構造を抽出します。例えば、ある企業の顧客データを教師なし学習で解析することで、似たような傾向を持つグループや、注目すべき特徴を抽出することができます。

一方、教師あり学習は、データにラベルや目的変数がある場合に用いられます。目的変数と予測値との誤差を最小化するようにモデルを構築し、未知のデータに対して予測を行います。例えば、画像認識や音声認識などが教師あり学習の例です。

カーネル密度推定は、教師なし学習の手法の一つで、データがどのような分布に従っているかを推定する方法です。具体的には、データ点の周りにカーネルと呼ばれる関数を置き、そのカーネルの大きさや位置を調整することで、データ点の密度を推定します。これにより、データの分布を把握することができます。

教師なし学習とは、機械学習の手法の一つである。「出力すべきもの」があらかじめ決まっていないという点で教師あり学習とは大きく異なる。データの背後に存在する本質的な構造を抽出するために用いられる。

教師あり学習は、その「出力すべきもの」も入力として与える手法であり、データの背後に存在する本質的な構造を抽出するよりむしろ、思い通りの出力を再現する機械の構成に用いられる。

具体的な例として以下のようなものがある。

・クラスター分析

・主成分分析

・ベクトル量子化

・自己組織化マップ

教師なし学習

カーネル密度推定

カーネル密度推定は、確率変数の確率密度関数を推定する方法の一つで、パラメータが少ないモデルや確率分布を使わない、ノンパラメトリックな手法です。

パラメータが少ないモデルや確率分布を使わないということは、カーネル密度推定がモデルフリー(model-free)な手法であることを示しています。つまり、事前に確率分布の形を決めておく必要がなく、データ自体から確率密度関数を推定することができます。

これに対して、パラメトリックな手法は、データを表す確率分布の形をあらかじめ決めておく必要があります。例えば、正規分布、ポアソン分布、指数分布などのように、あらかじめ確率分布の形が決められているモデルを使います。そして、そのモデルに含まれるパラメータ(平均値、分散、密度関数の形状のパラメータなど)を推定することで、確率密度関数を推定します。このような手法をパラメトリックな手法と呼びます。

ノンパラメトリックな手法は、データの形状について事前に仮定を置かないため、より柔軟なモデルを構築することができます。一方で、データが大量にある場合、計算量が非常に大きくなる可能性があるため、パラメトリックな手法よりも時間がかかることがあります。また、推定された確率密度関数の精度についても、パラメトリックな手法と比較して評価が難しいことがあります。

例えば、100m走のタイムは連続的な値をとる確率変数です。このとき、確率密度関数を求めることで、ある値をとる確率を知ることができます。しかし、確率密度関数は通常未知であるため、カーネル密度推定を使って推定する必要があります。

カーネル密度推定では、ある値の周りにカーネルと呼ばれる関数を置き、そのカーネルの形や大きさを調整して、周りにあるデータ点の密度を推定します。例えば、100m走のタイムが10秒台である選手が多い場合、その周りにカーネルを置くことで、10秒台でのタイムをとる確率が高いことを推定できます。

カーネル密度推定は、パラメトリックな手法よりも柔軟性があり、より現実的な分布を表現できる場合があるため、様々な分野で使われています。

カーネル密度推定は、統計学において、確率変数の確率密度関数を推定するノンパラメトリック手法のひとつ。エマニュエル・パルツェンの名をとってパルツェン窓とも。

カーネル密度推定

カーネル密度推定は、ある母集団から抽出された標本データが与えられた場合に、そのデータを元にして母集団全体の分布を推定する方法です。具体的には、標本データを元にカーネル関数を適用し、その関数を用いてデータのない範囲(外側)の値を推定します。たとえば、ある地域の住民の身長を測定し、そのデータを元にして、その地域の住民全体の身長分布を推定することができます。その際に、カーネル密度推定を使用して、身長が測定されていない人たちの身長を推定することができます。

大まかに言えば、ある母集団の標本のデータが与えられたとき、カーネル密度推定を使えばその母集団のデータを外挿できる。

ヒストグラムは、一様なカーネル関数によるカーネル密度推定量と見ることもできる。

カーネル密度推定

「カーネル」とは、統計学においていくつかの異なる意味で使われる言葉です。ノンパラメトリックな推定法において、カーネルは重み付け関数として使われます。つまり、データに対してカーネル関数を適用することで、近傍のデータに重みをつけて平滑化を行うことができます。

ノンパラメトリックな推定を行う場合、通常は使うカーネル関数に加えて、カーネルの幅を指定する必要があります。幅を大きくすると、平滑化が強くなり、カーネルの幅を小さくすると、データに近い点の重みが強くなります。

カーネル関数には、一様、三角、Epanechnikov(エパネチコフ)、quartic (biweight)、triweight、tricube、ガウシアン(ガウス関数)、コサインなどがあります。例えば、ガウス関数を用いると、平滑化したときの曲線がなめらかになるため、よく使われます。

ノンパラメトリック手法において、カーネルとは、ノンパラメトリックな推定手法に用いられる重み付け関数のことである。カーネルは、確率変数の確率密度関数を推定するためのカーネル密度推定や、確率変数の条件付き期待値を推定するカーネル回帰に用いられる。カーネルは時系列分析においては窓関数という名称で、ピリオドグラムによってスペクトル密度を推定するのに用いられる。その他の利用法としては、点過程の時間可変な強度の推定にも用いられる。そこでは窓関数(カーネル)は、時系列データとともに畳み込まれる。

ノンパラメトリックな推定を実行する際はふつう、(カーネル関数に加えて)カーネルの幅も指定されなければならない。

カーネル(統計学)
カーネル関数グラフ
一様{\displaystyle K(u)={\frac {1}{2}}\,\mathbf {1} _{\{|u|\leq 1\}}}Kernel uniform.svg
三角{\displaystyle K(u)=(1-|u|)\,\mathbf {1} _{\{|u|\leq 1\}}}Kernel triangle.svg
Epanechnikov
(エパネチコフ)
{\displaystyle K(u)={\frac {3}{4}}(1-u^{2})\,\mathbf {1} _{\{|u|\leq 1\}}}Kernel epanechnikov.svg
quartic
(biweight)
{\displaystyle K(u)={\frac {15}{16}}(1-u^{2})^{2}\,\mathbf {1} _{\{|u|\leq 1\}}}Kernel quartic.svg
triweight{\displaystyle K(u)={\frac {35}{32}}(1-u^{2})^{3}\,\mathbf {1} _{\{|u|\leq 1\}}}Kernel triweight.svg
tricube{\displaystyle K(u)={\frac {70}{81}}(1-{\left|u\right|}^{3})^{3}\,\mathbf {1} _{\{|u|\leq 1\}}}Kernel tricube.svg
ガウシアン
(ガウス関数)
{\displaystyle K(u)={\frac {1}{\sqrt {2\pi }}}\exp \left(-{\frac {1}{2}}u^{2}\right)}Kernel exponential.svg
コサイン{\displaystyle K(u)={\frac {\pi }{4}}\cos \left({\frac {\pi }{2}}u\right)\mathbf {1} _{\{|u|\leq 1\}}}Kernel cosine.svg
ロジスティック{\displaystyle K(u)={\frac {1}{e^{u}+2+e^{-u}}}}Kernel logistic.svg
Silverman
カーネル
{\displaystyle K(u)={\frac {1}{2}}\exp \left(-{\frac {|u|}{\sqrt {2}}}\right)\sin \left({\frac {|u|}{\sqrt {2}}}+{\frac {\pi }{4}}\right)}Kernel Silverman.svg
カーネル関数

定義

「確率密度変数ƒの標本」とは、ƒと同じ分布を持つ複数のデータの集まりのことです。 カーネル密度推定量は、このデータを元に、Kとhを使ってƒの推定を行う手法です。Kは重み付け関数で、データが密集している領域で高い値を持ち、データが疎である領域で低い値を持ちます。hは平滑化の度合いを調整するパラメータで、hが小さいときにはより詳細な解像度でƒを推定することができますが、ノイズが大きくなる可能性があります。逆に、hが大きすぎるとƒが滑らかになりすぎて、データの特徴を見失うことがあります。

{\displaystyle {\hat {f}}_{h}(x)={\frac {1}{nh}}\sum _{i=1}^{n}K\left({\frac {x-x_{i}}{h}}\right)}

カーネル関数としてよく使われるのは、平均が0で分散が1の標準ガウス関数です。

{\displaystyle K(x)={1 \over {\sqrt {2\pi }}}\,e^{-x^{2}/2}}

Pythonプログラミング

「カーネル密度推定」をイメージしやすいようPythonでのプログラミングについても学びます。

scikit-learn トイデータセット

機械学習ライブラリscikit-learnに用意されている「トイデータセット」を使います。トイデータセットは、機械学習の問題を解くためのサンプルデータセットのことで、いくつかの種類が用意されています。例えば、Iris(アヤメ)の花の特徴から、その種類を分類する問題を解くための「irisデータセット」や、ボストン市の住宅価格に関するデータを用いて、住宅価格を予測する問題を解くための「bostonデータセット」などがあります。

Irisデータセットは、Setosa、Versicolour、Virginicaの3種類のアヤメの花のデータが含まれており、各種類について50個のサンプルがあります。このデータセットには、がく片の長さや幅、花びらの長さや幅、そしてアヤメの種類という4つの特徴量が含まれています。

アヤメの花の特徴量(花弁の長さや幅、がく片の長さや幅)を読み込み、カーネル密度推定を行い、グラフにプロットすることができます。

scikit-learnにはいくつかの小さな標準データセットが付属しており、外部のウェブサイトからファイルをダウンロードする必要はありません。

これらは以下の関数を使って読み込むことができます。

load_boston() : load_boston は 1.0 で非推奨となり、1.2 で削除される予定である。

load_iris() : アヤメのデータセット(分類)をロードして返す。

load_diabetes() : 糖尿病のデータセット(回帰)をロードして返す。

load_digits() : 数字のデータセット(分類)をロードして返す。

load_linnerud() : 身体運動のデータセットをロードして返す。

load_wine() : ワインのデータセット(分類)をロードして返す。

load_breast_cancer() : ウィスコンシン州の乳がんのデータセット(分類)をロードして返す。7.1. Toy datasetsts

アヤメのデータセット

アヤメのデータセットは、機械学習における代表的なデータセットの1つであり、サンプルデータの中でも特に有名なものの1つです。

アヤメのデータセットには、ヒオウギアヤメ(Iris-Setosa)、アイリス・バージカラー(Iris Versicolour)、アイリス・ヴァージニカ(Iris Virginica)の3種類のアヤメが含まれており、各々50件のデータがあります。データには、「がく片の長さ」、「がく片の幅」、「花びらの長さ」、「花びらの幅」といった4つの属性があり、また3種のアヤメの分類に関するデータも含まれています。

Iris Setosa (ヒオウギアヤメ)

アラスカ、メイン、カナダ(ブリティッシュコロンビア、ニューファンドランド、ケベック、ユーコンなど)、ロシア(シベリアなど)、アジア北東部、中国、韓国、日本など北極海を越えて広く分布する根生葉の多年草です。

茎は高く伸び、葉は中緑色、花は紫、紫紺、青、ラベンダー色です。また、ピンクや白の花を咲かせる植物もあります。

Irissetosa1.jpg
Iris Setosa (ヒオウギアヤメ)

Iris Versicolour (アイリス・バージカラー)

北アメリカ、アメリカ東部とカナダ東部に自生するアヤメの一種です。スゲ草地や湿地、川岸や海岸に普通に見られます。

特異形質バージカラー(versicolor)は「様々な色彩の」という意味です。

Blue Flag, Ottawa.jpg
Iris Versicolour (アイリス・バージカラー)

Iris Virginica (アイリス・ヴァージニカ)

北アメリカ東部原産の多年草です。アメリカ南東部のフロリダ州からジョージア州にかけての海岸平野によく見られます。

Iris virginica 2.jpg
Iris Virginica (アイリス・ヴァージニカ)

プログラム

アヤメのデータセットのSepal lengthについて、3種類のアヤメのそれぞれについてカーネル密度推定を行い、グラフにプロットします。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from scipy.stats import gaussian_kde

# アヤメのデータセットを読み込む
iris = load_iris()

# 特徴量データをDataFrameに変換する
df = pd.DataFrame(iris.data, columns=iris.feature_names)

# Setosa、Versicolour、Virginicaの各データのインデックスを取得する
setosa_idx = iris.target == 0
versicolor_idx = iris.target == 1
virginica_idx = iris.target == 2

# 各データのSepal lengthのカーネル密度推定を計算する
setosa_density = gaussian_kde(df.loc[setosa_idx, 'sepal length (cm)'])
versicolor_density = gaussian_kde(df.loc[versicolor_idx, 'sepal length (cm)'])
virginica_density = gaussian_kde(df.loc[virginica_idx, 'sepal length (cm)'])

# x軸の範囲を決定する
x_min = df['sepal length (cm)'].min()
x_max = df['sepal length (cm)'].max()
x_range = np.linspace(x_min, x_max, 100)

# グラフを描画する
plt.plot(x_range, setosa_density(x_range), label='Setosa')
plt.plot(x_range, versicolor_density(x_range), label='Versicolor')
plt.plot(x_range, virginica_density(x_range), label='Virginica')
plt.legend()
plt.xlabel('Sepal length (cm)')
plt.ylabel('Density')
plt.show()

実行結果

以下のようなグラフが表示されます。

このグラフから、3種類のアヤメのうち、Setosaは他の2種類よりも短いSepal lengthを持つ傾向があることがわかります。一方、VersicolourとVirginicaは比較的似た範囲のSepal lengthを持っていますが、Virginicaの方が少し長い傾向があります。これは、カーネル密度推定を使用して、アヤメのデータセットから洞察を得る例となっています。

プログラムの説明

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from scipy.stats import gaussian_kde

必要なライブラリをインポートしています。NumPyは数値計算、Pandasはデータ処理、Matplotlibはグラフ描画、Scikit-learnのdatasetsからirisデータセットを読み込み、Scipyのstatsからカーネル密度推定のためのgaussian_kdeをインポートしています。

# アヤメのデータセットを読み込む
iris = load_iris()

# 特徴量データをDataFrameに変換する
df = pd.DataFrame(iris.data, columns=iris.feature_names)

Scikit-learnからirisデータセットを読み込んでいます。irisデータセットは、アヤメの花の種類(setosa, versicolor, virginica)ごとに、がく片の長さ(sepal length)、がく片の幅(sepal width)、花びらの長さ(petal length)、花びらの幅(petal width)の4つの特徴量を持った150のサンプルが含まれています。

irisデータセットの特徴量データをDataFrameに変換しています。DataFrameは表形式のデータ構造で、行と列のラベルがついています。

# Setosa、Versicolour、Virginicaの各データのインデックスを取得する
setosa_idx = iris.target == 0
versicolor_idx = iris.target == 1
virginica_idx = iris.target == 2

Setosa、Versicolour、Virginicaの各データのインデックス(行に対して割り当てられる一意のラベル)を取得しています。

アヤメの種類は、setosa、versicolor、virginicaの3つがあり、それぞれのインデックスを取得するためには、iris.targetの値が0、1、2である行をフィルタリングする必要があります。このプログラムでは、比較演算子「==」を使用して、Irisデータセットのtarget変数が0、1、2である行のインデックスを取得しています。

target属性を使用することで、各サンプルが属する花の種類を表す0, 1, 2のラベルが格納された配列を取得し、それぞれの種類に属するサンプルのインデックスを取得しています。

# 各データのSepal lengthのカーネル密度推定を計算する
setosa_density = gaussian_kde(df.loc[setosa_idx, 'sepal length (cm)'])
versicolor_density = gaussian_kde(df.loc[versicolor_idx, 'sepal length (cm)'])
virginica_density = gaussian_kde(df.loc[virginica_idx, 'sepal length (cm)'])

各種類のデータのSepal lengthのカーネル密度推定を計算しています。df.loc[setosa_idx, 'sepal length (cm)’]のように、locメソッドを使用してDataFrameから特定の行と列を抽出し、gaussian_kde関数を用いてカーネル密度推定を計算しています。

gaussian_kde()は、カーネル密度推定 (Kernel Density Estimation, KDE) を行うための関数です。KDEは、データから確率密度関数 (Probability Density Function, PDF) を推定する手法であり、その推定にカーネル関数と呼ばれる関数を用います。

具体的には、KDEは、各データ点の周囲にカーネル関数を重ね合わせたものを正規化することで、PDFを推定します。この際、カーネル関数は、例えば正規分布や三角形のような形状をしていることがあります。

gaussian_kde()は、KDEを計算する際に用いるカーネル関数として、ガウス関数を用います。ガウス関数は、正規分布の確率密度関数として知られており、以下の式で表されます。

ここで、xは入力値、μは平均値、σは標準偏差を表します。

gaussian_kde()は、引数に与えられたデータから、上記の式を用いてガウスカーネルを生成し、それらを重ね合わせてKDEを計算します。計算結果は、関数の形で表され、入力値に対する確率密度を出力します。

# x軸の範囲を決定する
x_min = df['sepal length (cm)'].min()
x_max = df['sepal length (cm)'].max()
x_range = np.linspace(x_min, x_max, 100)

x軸の範囲を決定しています。x_minとx_maxは、DataFrameの’sepal length (cm)’列の最小値と最大値を取得しています。x_rangeは、x_minからx_maxまで100点に分割した等差数列を生成しています。

# グラフを描画する
plt.plot(x_range, setosa_density(x_range), label='Setosa')
plt.plot(x_range, versicolor_density(x_range), label='Versicolor')
plt.plot(x_range, virginica_density(x_range), label='Virginica')
plt.legend()
plt.xlabel('Sepal length (cm)')
plt.ylabel('Density')
plt.show()

3つの異なるアヤメの品種(Setosa、Versicolour、Virginica)について、Sepal lengthのカーネル密度推定を計算し、それらを1つのグラフにプロットしています。

具体的には、NumPyのlinspace()関数を使用して、x_range変数にdf['sepal length (cm)’]列の最小値と最大値を指定した100個の等間隔の数値を作成します。その後、setosa_density()、versicolor_density()、およびvirginica_density()の3つのカーネル密度推定関数を使用して、それぞれの品種についてSepal lengthの密度分布を計算します。

最後に、Matplotlibのplot()関数を使用して、x軸にはx_rangeを、y軸には各品種の密度分布をプロットします。legend()関数で品種を凡例に追加し、xlabel()関数とylabel()関数で軸ラベルを設定します。そして、show()関数でグラフを表示します。