ノンパラメトリック密度推定

ノンパラメトリック密度推定#

2つの正規分布からなる混合分布があるとする

0.3N(0.25,0.01)+0.7N(0.75,0.01)

確率密度関数は下の図の左側のようになる。

この分布から100個のサンプルが得られたとする。ヒストグラムは右側の図のようになる。

Hide code cell source 
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

# gen data
# 2つの正規分布の混合分布
w = 0.3
params = [dict(loc=0.25, scale=0.1), dict(loc=0.75, scale=0.1)]

def pdf(x):
    return w * norm.pdf(x, **params[0]) + (1 - w) * norm.pdf(x, **params[1])


def rvs(size):
    n = int(w * size)
    return np.append(norm.rvs(size=n, random_state=0, **params[0]), norm.rvs(size=(size - n), random_state=0, **params[1]))


fig, axes = plt.subplots(figsize=[10, 3], ncols=2)
axes[0].set(xlabel="x", ylabel="probability density", title="probability density function")
x_range = np.linspace(0, 1, 100)
axes[0].plot(x_range, [pdf(xi) for xi in x_range], color="dimgray")

axes[1].set(xlabel="x", ylabel="count", title="histogram")
n = 50
x = rvs(n)
x = np.sort(x)
axes[1].hist(x, color="dimgray")
fig.show()
../_images/84ed31d6f9bb73fa850a9103215914ae792d62b70a42b63cb5a6b0cc322c6c61.png

ヒストグラム密度推定法#

ヒストグラム密度推定法（histogram density estimation method）

1つの連続変数xが対象の場合を考える。

標準的なヒストグラムでは、xを幅Δiの区間に区切り、i番目の区間に入ったxの観測値の数niと観測値の総数Nを用いて、各区間の確率密度を

pi=niNΔi

で推定する。

Hide code cell source 
class HistogramEstimator:

    def __init__(self, bandwidth: float):
        self.h = bandwidth

    def fit(self, x: np.array):
        # calc bins
        range_ = max(x) - min(x)
        times = int(np.ceil(range_ / self.h))
        bins = [min(x) + self.h * i for i in range(times + 1)]
        # calc probabilities
        N = x.shape[0]
        probs = []
        deltas = []
        for i in range(1, len(bins)):
            values = x[(bins[i - 1] <= x) & (x < bins[i])]
            n = len(values)
            delta = bins[i] - bins[i - 1]
            p = n / (N * delta)
            probs.append(p)
            deltas.append((bins[i - 1], bins[i]))
        self.probs_ = probs
        self.deltas_ = deltas
        return self

    def score(self, x: np.array) -> np.array:
        p_pred = np.zeros(shape=x.shape[0])
        for (lower, upper), p in zip(self.deltas_, self.probs_):
            p_pred[(lower <= x) & (x < upper)] = p
        return p_pred

# plot
fig, axes = plt.subplots(figsize=[14, 3], ncols=3)

for ax, h in zip(axes, [0.05, 0.10, 0.2]):
    # true
    x_range = np.linspace(0, 1, 100)
    ax.plot(x_range, [pdf(xi) for xi in x_range], color="dimgray", label="True PDF")

    # estimates
    estimator = HistogramEstimator(bandwidth=h)
    estimator.fit(x)
    for (lower, upper), p in zip(estimator.deltas_, estimator.probs_):
        ax.bar(lower, p, h, color="white", edgecolor="steelblue")
        
    ax.bar(lower, p, h, color="white", edgecolor="steelblue", label="histogram method")

    ax.scatter(x, np.zeros(x.shape[0]), label="samples", color="gray", alpha=0.5)
    ax.set(xlabel="x", ylabel="probability density", title=f"Histogram method (h={h})")
    ax.legend()

fig.show()
../_images/a3e18a857bd48b4565858c0638f02fed3b3f73d47cee3ed834040e9289c88a92.png

この方法の結果の良し悪しは区間幅Δiに大きく依存する。幅が狭すぎても区間に含まれるサンプルが少なすぎて推定のばらつきが大きくなるし、幅が広すぎても表現力が不足して分布をうまく捉えられなくなる。

また、ヒストグラム法の問題として

  1. 推定した密度が区間の縁で不連続になる

  2. 次元の呪いに弱い：次元数を上げていった場合、D次元空間を各変数につきM個の区間にすると、区間の総数はMD個になり、各区間に含まれるデータ量が不足する

といったものがある

カーネル密度推定法#

D次元のユークリッド空間中の未知の確率密度p(x)から観測値の集合が得られていて、この集合からp(x)の値を推定したいとする。

xを含むある小さな領域Rを考える。この領域に割り当てられた確率は

P=∫Rp(x)dx

と表すことができる。

（参考）領域Rのイメージ（1次元の場合）
../_images/5f5862170226b62b3bc9569fbe0f2156288108575c34c4faacfa0e82fd6ee570.png

ここでp(x)から得られたN個の観測値からなるデータ集合があるとする。各データ店が領域R中にある確率はPなので、R内の点の総数Kは二項分布に従う

Bin(K|N,P)=N!K!(N−K)!PK(1−P)N−K

よって、データ点がこの領域内にある平均割合と分散は

E[K/N]=P,Var[K/N]=P(1−P)N

となる。

大きいNについては、分散が小さくなって平均の周囲で鋭く尖った分布となり、

K≃NP

となる。

Rが確率密度p(x)がこの領域内でほぼ一定とみなせるほど十分に小さいものであると仮定できるのであれば、領域の体積Vを用いて

P≃p(x)V

となる。

これらを組み合わせて、次の密度の推定量が得られる。

p(x)=KNV

Note

上記の推定量の正しさは、2つの相反する仮定に基づく

  1. 領域内では密度が一定とみなせるほど十分に領域Rは小さい

  2. 領域内のデータ点Kが二項分布が尖るほどに多く存在している

カーネル関数#

確率密度を求めたいデータ点xを中心とする小さな超立方体を領域Rとする。 この領域内にある点の数Kを数えるには、次の関数を定義しておくと便利である。

k(u)={1,|ui|≤12, if i=1,⋯,D0,otherwise

これは原点を中心とする単位立方体を表す。 関数k(u)はカーネル関数（kernel function）のひとつであり、今回の用途ではParzen窓（parzen window）とも呼ばれる。

Hide code cell source 
def k(u):
    return 1 * all([abs(ui) <= (1/2) for ui in u])

u1 = np.linspace(-1, 1, 51)
u2 = np.linspace(-1, 1, 51)
u1, u2 = np.meshgrid(u1, u2)
k_ = np.array([k([_u1, _u2]) for _u1, _u2 in zip(u1.flatten(), u2.flatten())]).reshape(u1.shape)

fig, ax = plt.subplots(dpi=100, subplot_kw={"projection": "3d"})
ax.plot_surface(u1, u2, k_, rstride=1, cstride=1, linewidth=1, antialiased=True, alpha=0.7)
ax.set(xlabel="u1", ylabel="u2", zlabel="k", title="k(u)", zlim=[0, 2])
fig.show()
../_images/07927ed45e45e8f8e7b65650dda54fa1a9f7ae3845508fa159ad5114773856a0.png

k((x−xn)/h)はxを中心とする一辺がhの立方体の内部に、データ点xnがあれば1に、そうでなければ0となる。

例えばx=(2,3),h=2の場合は次の図のようになる

Hide code cell source 
_x = np.array([2, 3])
h = 2

x1 = np.linspace(-1, 5, 51)
x2 = np.linspace(0, 5, 51)
x1, x2 = np.meshgrid(x1, x2)
k_ = np.array([k( ( _x - np.array([_x1, _x2]) ) / h )
               for _x1, _x2 in zip(x1.flatten(), x2.flatten())]).reshape(u1.shape)

fig, ax = plt.subplots(dpi=100, subplot_kw={"projection": "3d"})

ax.scatter(_x[0], _x[1], color="orange")
ax.text(_x[0] * 1.1, _x[1] * 1.1, 0, f"x=({_x[0]}, {_x[1]})", color="orange")

ax.plot_surface(x1, x2, k_, rstride=1, cstride=1, linewidth=1, antialiased=True, alpha=0.4)
ax.set(xlabel="x1", ylabel="x2", zlabel="k", title="k([x - x_n] / h)", zlim=[0, 2.5])
fig.show()
../_images/c1506f121bbf6e85111a713cddffab778733b62a0bad839368b5309b487ba4a2.png

この立方体内部の総点数は

K=∑n=1Nk(x−xnh)

となる。

さきほどのp(x)の推定量

p(x)=KNV

に代入すると

p(x)=1NV∑n=1Nk(x−xnh)

一辺がhのD次元立方体の体積がV=hDであることを用いると

p(x)=1N1hD∑n=1Nk(x−xnh)

となる。

このカーネルを使用した推定結果は次の図のようになる。 立方体を重ねるような推定を行うため平滑性がなく、ギザギザした密度関数が推定されている。

Hide code cell source 
class KernelEstimator:

    def __init__(self, bandwidth: float):
        self.h = bandwidth

    def fit(self, X: np.array):
        self.X_train_ = X
        self.N_ = X.shape[0]
        self.D_ = X.shape[1]
        return self

    def score(self, X: np.array) -> np.array:
        p_preds = []
        for x in X:
            U = (x - self.X_train_) / self.h
            K = [k(u) for u in U]
            p_pred = (1 / self.N_) * (1 / self.h ** self.D_) * np.sum( K )
            p_preds.append(p_pred)
        return np.array(p_preds)

# plot
fig, axes = plt.subplots(figsize=[14, 3], ncols=3)

for ax, h in zip(axes, [0.05, 0.10, 0.20]):
    # true
    x_range = np.linspace(0, 1, 100)
    ax.plot(x_range, [pdf(xi) for xi in x_range], color="dimgray", label="True PDF")

    # estimates
    X_train = x.reshape(-1, 1)
    X_test = x_range.reshape(-1, 1)

    estimator = KernelEstimator(bandwidth=h)
    estimator.fit(X_train)
    p = estimator.score(X_test)
    ax.plot(x_range, p, label="estimated PDF")

    ax.scatter(x, np.zeros(x.shape[0]), label="samples", color="gray", alpha=0.5)
    ax.set(xlabel="x", ylabel="probability density", title=f"Kernel method (h={h})")
    ax.legend()

fig.show()
../_images/09c9b5abf533d3a526125a6fad654ab8bcb3926bfbb970ca8466f2ed92ad40f1.png

ガウシアンカーネル#

ガウス分布（正規分布）をカーネル関数に用いることで滑らかな密度推定を行う。

p(x)=1N∑n=1N1(2πh2)D/2exp⁡{−||x−xn||22h2}
Hide code cell source 
from scipy.stats import norm

class GaussianKernelEstimator:

    def __init__(self, bandwidth: float):
        self.h = bandwidth

    def fit(self, X: np.array):
        self.X_train_ = X
        self.N_ = X.shape[0]
        self.D_ = X.shape[1]
        return self

    def score(self, X: np.array) -> np.array:
        p_preds = []
        for x in X:
            p_pred = norm.pdf(x=x, loc=self.X_train_, scale=self.h).mean()
            p_preds.append(p_pred)
        return np.array(p_preds)

# plot
fig, axes = plt.subplots(figsize=[14, 3], ncols=3)

for ax, h in zip(axes, [0.02, 0.05, 0.10]):
    # true
    x_range = np.linspace(0, 1, 100)
    ax.plot(x_range, [pdf(xi) for xi in x_range], color="dimgray", label="True PDF")

    # estimates
    X_train = x.reshape(-1, 1)
    X_test = x_range.reshape(-1, 1)

    estimator = GaussianKernelEstimator(bandwidth=h)
    estimator.fit(X_train)
    p = estimator.score(X_test)
    ax.plot(x_range, p, label="estimated PDF")

    ax.scatter(x, np.zeros(x.shape[0]), label="samples", color="gray", alpha=0.5)
    ax.set(xlabel="x", ylabel="probability density", title=f"Gaussian Kernel method (h={h})")
    ax.legend()

fig.show()
../_images/d3d337e3ab955869949a48370dc5177003fef9406741996a12d54a2f68f64710.png
Hide code cell source 
# 数式からガウスカーネルを実装したけど合わなかったやつ
def gauss(x, x_n, h, D):
    y = (1 / (2 * np.pi * h **2)**(D/2)) * np.exp(-(np.linalg.norm(x - x_n)**2 / 2 * h ** 2))
    return np.mean(y)

class GaussianKernelEstimator:

    def __init__(self, bandwidth: float):
        self.h = bandwidth

    def fit(self, X: np.array):
        self.X_train_ = X
        self.N_ = X.shape[0]
        self.D_ = X.shape[1]
        return self

    def score(self, X: np.array) -> np.array:
        p_preds = []
        for x in X:
            p_pred = gauss(x, x_n=self.X_train_, h=self.h, D=self.D_)
            p_preds.append(p_pred)
        return np.array(p_preds)

# plot
fig, axes = plt.subplots(figsize=[14, 3], ncols=3)

for ax, h in zip(axes, [0.05, 0.10, 0.20]):
    # true
    x_range = np.linspace(0, 1, 100)
    ax.plot(x_range, [pdf(xi) for xi in x_range], color="dimgray", label="True PDF")

    # estimates
    X_train = x.reshape(-1, 1)
    X_test = x_range.reshape(-1, 1)

    estimator = GaussianKernelEstimator(bandwidth=h)
    estimator.fit(X_train)
    p = estimator.score(X_test)
    ax.plot(x_range, p, label="estimated PDF")

    ax.scatter(x, np.zeros(x.shape[0]), label="samples", color="gray", alpha=0.5)
    ax.set(xlabel="x", ylabel="probability density", title=f"Kernel method (h={h})")
    ax.legend()

fig.show()
Hide code cell output
../_images/45d4cb3b6ee7103094e2878953b76aa179293dd2067818ea382b19bfade2c3bb.png

（参考）Scikit-learn実装#

2.8. Density Estimation — scikit-learn 1.2.2 documentation

（※kernel='tophat'は立方体カーネルに近いk(u;h)∝1 if u<h というもの）

Hide code cell source 
from sklearn.neighbors import KernelDensity

fig, axes = plt.subplots(figsize=[14, 3], ncols=3)
for ax, h in zip(axes, [0.02, 0.05, 0.10]):
    # true
    x_range = np.linspace(0, 1, 100)
    ax.plot(x_range, [pdf(xi) for xi in x_range], color="dimgray", label="True PDF")
    ax.scatter(x, np.zeros(x.shape[0]), label="samples", color="gray", alpha=0.5)

    # estimates
    X_train = x.reshape(-1, 1)
    X_test = x_range.reshape(-1, 1)
    
    for method in ["gaussian", "tophat"]:
        kde = KernelDensity(kernel=method, bandwidth=h).fit(X_train)
        log_likelihood = kde.score_samples(X_test)
        ax.plot(x_range, np.exp(log_likelihood), label=f"kernel='{method}'")
    ax.set(xlabel="x", ylabel="probability density", title=f"Kernel method (h={h})")
    ax.legend()
fig.show()
../_images/9a5f392038038ccffdbd81e41bb635050fd2d550c4c2fabab8e1de8610872eeb.png
Contents