ガウス過程を用いた回帰

概要

こんにちは、kwashiです。昨今、機械学習の流行りとともに、回帰モデルを用いた予測を制御に利用するなど、産業分野にも適用されるようになってきています。しかし、製品が関わる重要な場面では、予測の不確かさも重要になってきます。例えば、工場の環境(温度など）から、生産装置のパラメータを回帰するとします。しかし、学習に使用していないようなパターンが不意に発生しても、何らかのパラメータを回帰してしまいます。そこで、回帰はどの程度信頼できるのか、推論の不確かさも重要になってきます。

この不確かさを確率分布で示す分野にベイズ統計があります。ベイズ統計は、自動運転や音響処理など様々な分野に使われています。本記事では、ベイズ統計を用いた回帰の例を示します。ターゲットとして単純なモデルであるsin関数を回帰させます。最初に多項式回帰の例を示します。しかし、多項式では、パラメータを十分に設定しなければsin関数を表現できませんでした。そこで、ガウス過程を用いて関数そのものを推定する手法に発展させます。

また、本記事では、pythonのPyMCとScipyライブラリを中心に作成しています。

ベイズ統計に関して、本記事では十分にカバーできていません。私は以下の書籍を参考に本記事を書いたので、ぜひ読んでみてください。
続・わかりやすいパターン認識―教師なし学習入門
Bayesian Analysis with Python

Sin関数の多項式回帰

本章では、一様分布からランダムに抽出した角度xにおいて、平均sin(x) 標準偏差0.2である正規分布に則ってデータを生成する。また、最小二乗法を用いて、5次多項式にフィッテングさせる。5次多項式は、以下の式です。

$$f(x) = b_0 + b_1 x + b_2 x^2 + b_3 x^3 + b_4 x^4 + b_5 x^5$$

まず、結果を以下に示す。緑線がsin関数で、青点が学習に用いるデータである。また、赤が、単純に5次多項式を最小二乗法でフィッテングさせた結果である。

プログラムは以下のとおりです。scipyを用いて最小二乗法によるフィッテングを行っています。

import numpy as np
from scipy.odr import odrpack as odr
from scipy.odr import models
import matplotlib.pyplot as plt
import theano.tensor as tt
import scipy.stats as stat
import pymc3 as pm

# Sin 関数
np.random.seed(5)
x = np.random.uniform(0, 10, size=20)
y = np.random.normal(np.sin(x), 0.2) #正規分布, (平均, 分散)

#scipyによるフィッティング
polyFunc = models.polynomial(5) #5次
data = odr.Data(x, y)
myodr = odr.ODR(data, polyFunc, maxit=200)

myodr.set_job(fit_type=2) #最小二乗法

fit = myodr.run()

coeff = fit.beta[::-1]
err = fit.sd_beta[::-1]
print(coeff) #[-2.48756876e-05 -1.07573181e-02  2.09111045e-01 -1.21915241e+00, 2.11555200e+00 -1.08779827e-01] #x50

ベイズによる多項式回帰

本章では、前章で生成したデータをもとに、ベイズの枠組みで5次多項式を求めます。
モデルは、以下のプログラムのように、

• 正規分布より生成されたパラメータを持つ多項式から生成される平均
• コーシ分布から生成されるデータを分散
• 指数関数から生成されるデータを自由度 (正規分布の中央に集まる具合を操作するパラメータ)

とした、スチューデントのt分布

$$
{\rm St}(x|\mu,\lambda,\nu)={\Gamma(\nu/2+1/2)\over\Gamma(\nu/2)}\left({\lambda\over\pi\nu}\right)^{1/2}\left[1+{\lambda(x-\mu)^2\over\nu}\right]^{-\nu/2-1/2}
$$
とする。
そして、マルコフ連鎖モンテカルロ法(Markov chain Monte Carlo methods; MCMC)を用いて、パラメータを学習する。

with pm.Model() as model:
  b0 = pm.Normal('b0', mu=0, sd=10)
  b1 = pm.Normal('b1', mu=0, sd=5)
  b2 = pm.Normal('b2', mu=0, sd=5)
  b3 = pm.Normal('b3', mu=0, sd=5)
  b4 = pm.Normal('b4', mu=0, sd=5)
  b5 = pm.Normal('b5', mu=0, sd=5)

  #コーシー分布 x ∈ [0, ∞) https://docs.pymc.io/api/distributions/continuous.html#pymc3.distributions.continuous.HalfCauchy
  sd = pm.HalfCauchy('sd', 5)   
  mu = b0 + b1 * x + b2 * x**2 + b3 * x**3 + b4 * x**4 + b5 * x**5
  nu = pm.Exponential('nu', 1/30)

  y_pred = pm.StudentT('y_pred', mu=mu, sd=sd, nu=nu, observed=y)
  step = pm.Metropolis()
  trace = pm.sample(100000,step=step)

推定したパラメータを用いた多項式を赤で以下に示す。緑線がsin関数で、青点が学習に用いるデータ、青が最小二乗法の結果です。また、MCMCの結果得られた分散の平均が0.23で、生成モデルの分散が推定できていた。しかし、実際のsin関数と誤差が大きく、x=10付近などデータが少ない部分でも、分散が一定の0.23から生成される多項式になることが見て取れる。これは、モデルの設計が悪い。

ガウス過程への発展

前章では、関数f(x)を多項式のモデルで表現し、そのパラメータを推論していました。しかし、上手くデータを表現できておらず、モデルが間違っている可能性が大いにあります。そこで、関数を直接推論する方法を用います。
つまり、前章では、p(θ|x)のパラメータ推定を行っていましたが、本章では、p(f|x)を推論します。この関数の推論する手法としてガウス過程(GP)を利用します。

GPの公式な定義は、次の通りです。連続空間におけるすべての点は正規分布に従う変数と関連しており、GPはそれらの無限に多い確率変数の同時分布である。

しかし、ここでは、特に難しい議論をするつもりはありません。重要なことは、以下の3つです。

1. 事前分布が多変量正規分布である。f(x) ~ MvNormal(u|K)
2. 共分散行列(K)がカーネルである。
3. 共分散行列のパラメータ
4. 事後分布が、解析的に計算可能である。

2のカーネルは、サポートベクターマシンなどで用いられるカーネルを想像してください。カーネルは様々あり、共分散行列として、どのカーネルも使用可能です。本記事では、ガウシアンカーネルを用います。それでは、1, 2より得られるGPの事前分布を見てみましょう。以下のプログラムで、多変量正規分布(平均=0, 共分散行列=ガウシアンカーネル)をプロットしています。図を見ていると、様々な関数を表現していることがわかります。このような表現が無限にあります。
※注意点としては、別のカーネルを用いた場合、もとギザギザした関数など表現が変わってきます。しかし、今回は、sin関数を回帰すると考えたとき、以下の図のような表現が適切だと見て取れます。

def squareDistance(x, y):
  return np.array([[(x[i]-y[j])**2 for i in range(len(x))] for j in range(len(y))])

np.random.seed(1)
tp = np.linspace(0, 10, 100)

cov = np.exp(-squareDistance(tp, tp))
plt.plot(tp, stats.multivariate_normal.rvs(cov=cov, size=6).T

上図に示したような多変量正規分布は、共分散行列を変化させることで、様々な表現をしています。その共分散行列をベイズの枠組みで推定するためには、3で示したように推定すべきパラメータを設定する必要があります。そこで、共分散行列を

K_{i,j} = \eta \exp(-\rho D)  \ \  (i \neq j) \\
K_{i,j} = \eta + \sigma \ \ (i = j)

とします。

最後に、4解析的な事後分布の計算です。尤度を正規分布とすると、GPの事後分布は、解析的に問題を解くことが可能なため、GPの事後分布である多変量正規分布を単純な式で計算できます。事後分布により、回帰関数に予測が可能です。

詳しい式の説明は、[ガウス過程と機械学習]ガウス過程回帰に記載されています。

以下に、予測結果とプログラムの一部を示します。推定したGPの事後分布より予測した結果が赤線で、緑線がsin関数で、青点が学習に用いるデータ、青線が最小二乗法の結果です。GPで予測した結果、赤線が青線より緑線に近いことから、パラメータの調整を行わず、直接性能が良い回帰モデルが推定できていることがわかります。また、赤線の散らばり具合を見たところ、10付近などデータがないところでは、線が散らばる、つまり分散が大きいことが見て取れます。

#GP 事前分布: f(x) ~ GP(u=(u1, u2, ...), k(x, x')), u1, u2 .. = 0と仮定
#尤度: p(y|x, f(x)) ~ N(f, σ^2I)
#GP 事後分布 p(f(x)|x,y) ~GP(U,∑)

def squareDistance(x, y):
  return np.array([[(x[i]-y[j])**2 for i in range(len(x))] for j in range(len(y))])

points = np.linspace(xmin, xmax, xsize) #参照点

with pm.Model() as model:
  #事前分布
  mu = np.zeros(xsize)
  eta = pm.HalfCauchy('eta', 5)
  rho = pm.HalfCauchy('rho', 5)
  sigma = pm.HalfCauchy('sigma', 5)

  D = squareDistance(x, x)
  K = tt.fill_diagonal(eta * pm.math.exp(-rho * D), eta + sigma) #i=jの位置にeta+sigmaを格納

  gpPri = pm.MvNormal('obs', mu, cov=K, observed=y)

  #事後分布の計算(ガウス分布の事後分布は、式で計算可能)
  K_ss = eta * pm.math.exp(-rho * squareDistance(points, points))
  K_s = eta * pm.math.exp(-rho * squareDistance(x, points))

  MuPost = pm.Deterministic('muPost', pm.math.dot(pm.math.dot(K_s, tt.nlinalg.matrix_inverse(K)), y)) #K_s.T inv(K) y
  SigmaPost = pm.Deterministic('sigmaPost', K_ss - pm.math.dot(pm.math.dot(K_s, tt.nlinalg.matrix_inverse(K)), K_s.T))

  step = pm.Metropolis()
  trace = pm.sample(10000,step=step) 

まとめ

本記事では、不確かさの表現が重要であることからベイズ統計を導入しました。しかし、sin関数の回帰モデルを作成する上で、sin関数を上手く表現してくれるモデルを設計することが難しかったことから、ガウス過程を導入しました。詳しい説明を避け、概念を説明したつもりです。詳しい内容は、記事に載せている他記事や本を参考にしてください。