Kunst Site
Blogブログ
About Site当サイトについて
Portfolioポートフォリオ
Contact問い合わせ

(実装)ガウス過程回帰を理解したい-Gaussian Process Regression


データサイエンス教師あり学習
更新日:2022-02-22

本記事ではガウス過程回帰の実装を行っていきたいと思います。

理論などのインプットばかりではわかった気にはなっても使いこなすことはできないので、本記事を通して実装をアウトプットしました。

ガウス過程回帰とは

ガウス過程回帰は入力XXXから出力yyyへの関数を確率変数として推定します。ガウス過程はデータの平均値を予測すると同時に,その分散を出力します。

いま、学習データx1,⋯ ,xNx_1,\cdots,x_Nx1​,⋯,xN​に対応する出力f(x1),⋯ ,⋯ ,f(xN)f(x_1),\cdots,\cdots,f(x_N)f(x1​),⋯,⋯,f(xN​)がガウス分布に従っているとき、 新しい入力データ集合XXXに対応する出力yyyを予測したいときを考えると、その出力は次の多変量ガウス分布に従っています。

p(y∗∣x∗,X,Y)=N(k∗TK−1Y,k∗,∗−k∗TK−1k∗)p(y^{\ast}|\mathbf{x}^{\ast},X,Y)=\mathcal{N}(\mathbf{k}^T_{\ast}\mathbf{K}^{-1}\mathbf{Y},\mathbf{k}_{\ast,\ast}-\mathbf{k}^T_{\ast}\mathbf{K}^{-1}\mathbf{k}_{\ast})p(y∗∣x∗,X,Y)=N(k∗T​K−1Y,k∗,∗​−k∗T​K−1k∗​)

k∗\mathbf{k}_{\ast}k∗​は学習データXXXと新規のデータx∗\mathbf{x}^{\ast}x∗のカーネル、K\mathbf{K}Kは学習データXXXと学習XXXのカーネル、K∗,∗\mathbf{K}_{\ast,\ast}K∗,∗​は新規のデータx∗\mathbf{x}^{\ast}x∗と新規のデータx∗\mathbf{x}^{\ast}x∗のカーネルです.

使用データ

import numpy
np.random.seed(0)
resolution =  20
x_train = np.sort(np.random.rand(resolution) * 15)
y = np.sin(x_train) + np.random.rand(len(x_train))
x_test = np.linspace(0, 20, resolution)

xxxを0から150から150から15の範囲でランダムに生成してy=sin(x)+ϵy=sin(x)+\epsilony=sin(x)+ϵの関数へ入力して学習データを作成します。ϵ\epsilonϵはガウスノイズです。 テストデータに0から200から200から20の範囲を等間隔で学習データ分を入力として作成します。

学習データをプロットすると次のようになります。

学習データ

平均 k∗TK−1Y\mathbf{k}^T_{\ast}\mathbf{K}^{-1}\mathbf{Y}k∗T​K−1Y を求める

ガウス過程回帰では任意のカーネル関数を用います。予測される確率分布の平均k∗TK−1Y\mathbf{k}^T_{\ast}\mathbf{K}^{-1}\mathbf{Y}k∗T​K−1Yにもカーネル関数があります。

ここではカーネル関数はガウス関数を用います。

k(x,x′)=exp⁡(−σ2∣∣x−x′∣∣2)k(\mathbf{x},\mathbf{x'}) = \exp(-\frac{\sigma}{2}||\mathbf{x}-\mathbf{x'}||^2)k(x,x′)=exp(−2σ​∣∣x−x′∣∣2)
X_X_dist = np.sqrt((x_train[:, None] - x_train[None, :]) ** 2)
K = np.exp((-0.5 * self.σ * X_X_dist))
K_inv = np.linalg.inv(K)

X_test_dist = np.sqrt((x_train[:, None] - x_test[None, :]) ** 2)
k_star = np.exp((-0.5 * self.σ * X_test_dist))

y_pred_mean = k_star.T @ K_inv @ y

分散k∗,∗−k∗TK−1k∗\mathbf{k}_{\ast,\ast}-\mathbf{k}^T_{\ast}\mathbf{K}^{-1}\mathbf{k}_{\ast}k∗,∗​−k∗T​K−1k∗​を求める

予測される確率分布の分散k∗,∗−k∗TK−1k∗\mathbf{k}_{\ast,\ast}-\mathbf{k}^T_{\ast}\mathbf{K}^{-1}\mathbf{k}_{\ast}k∗,∗​−k∗T​K−1k∗​にもカーネル関数があります。

ここでのカーネル関数もガウス関数を用います。

k(x,x′)=exp⁡(−σ2∣∣x−x′∣∣2)k(\mathbf{x},\mathbf{x'}) = \exp(-\frac{\sigma}{2}||\mathbf{x}-\mathbf{x'}||^2)k(x,x′)=exp(−2σ​∣∣x−x′∣∣2)
X_X_dist = np.sqrt((x_train[:, None] - x_train[None, :]) ** 2)
K = np.exp((-0.5 * self.σ * X_X_dist))
K_inv = np.linalg.inv(K)

X_test_dist = np.sqrt((x_train[:, None] - x_test[None, :]) ** 2)
k_star = np.exp((-0.5 * self.σ * X_test_dist))
test_test_dist = np.sqrt((x_test - x_test) ** 2)
k_star_star = np.exp((-0.5 * self.σ * test_test_dist))

y_pred_cov = k_star_star - ((k_star.T @ K_inv) @ k_star)  # 分散共分散行列
y_pred_std = np.sqrt(np.diag(y_pred_cov))

ここで標準偏差=分散標準偏差=\sqrt{分散}標準偏差=分散​から標準偏差を求めています。

この標準偏差を求めることから、予測分布の不確かさがわかります。

実行結果

ガウス過程回帰の実装コードを描画のコードを足して Class にしてまとめました。

全体コードは次のようになります。

import numpy as np
import matplotlib.pyplot as plt

class gpr(object):
    def __init__(self, σ):
        self.σ = σ

    def fit(self, x_train, x_test, y):
        X_X_dist = np.sqrt((x_train[:, None] - x_train[None, :]) ** 2)
        K = np.exp((-0.5 * self.σ * X_X_dist))
        K_inv = np.linalg.inv(K)

        X_test_dist = np.sqrt((x_train[:, None] - x_test[None, :]) ** 2)
        k_star = np.exp((-0.5 * self.σ * X_test_dist))

        test_test_dist = np.sqrt((x_test - x_test) ** 2)
        k_star_star = np.exp((-0.5 * self.σ * test_test_dist))

        y_pred_mean = k_star.T @ K_inv @ y
        y_pred_cov = k_star_star - ((k_star.T @ K_inv) @ k_star)  # 分散共分散行列
        y_pred_std = np.sqrt(np.diag(y_pred_cov))  # 標準偏差
        return y_pred_mean, y_pred_cov, y_pred_std


if __name__ == "__main__":
    np.random.seed(0)
    resolution =  20
    x_train = np.sort(np.random.rand(resolution) * 15)
    y = np.sin(x_train) +  np.random.rand(len(x_train))

    x_test = np.linspace(0, 20, resolution)

    model = gpr(σ=1.0)
    mu, var, std = model.fit(x_train, x_test, y)

    # 描画
    fig = plt.figure(figsize=(12, 4))
    plt.cla()
    plt.scatter(x_train, y, c="black", marker="o", label="train data")
    plt.plot(x_test, mu, color="red", linewidth=1, label="mean")
    plt.fill_between(x_test, mu + std, mu - std, facecolor='green', alpha=0.2, label="confidence")
    plt.legend()
    plt.show()
結果

黒点は学習データです。赤い曲線が予測結果(推定したガウス分布の平均)で、緑色の領域は推定したガウス分布の共分散行列を使って計算した2σ2\sigma2σ区間を表していて、予測の信頼度が表現されています。

学習データの少ないところは分散(σ\sigmaσの値)が大きくなっているなど、関数の不確定性まで扱えています。

このように、予測値の不確実性を出力できることはガウス過程をはじめとするベイズモデルの長所です.

終わりに

ここまで読んでくださってありがとうございます。

編集リクエストもお待ちしてます。

Profile

kunst

データサイエンティスト兼フロントエンジニアのkunstです!
最近は山登りにハマっています。

目次

kunst Site

Welcome to kunst Site. Stay updated with the latest posts.

Privacy PolicyAbout Blog|BlogAbout SitePortfolioContact

Copyright © 2024 kunst Site