貝葉斯分析 高斯過程

2021-08-21 20:58:20 字數 3604 閱讀 6157

貝葉斯框架下, 可以用高斯過程來估計乙個函式 f : r→r. 對於每個xi, f(xi)可以用乙個均值方差暫未知的高斯分布來建模。因為連續空間的xi可以有無限個,擬合乙個函式的高斯過程其實乙個無限維的多元高斯。實際中,不管是我們的給定資料,還是測試點的個數都是有限的。因此無論是高斯過程先驗還是還是高斯過程後驗都是有限維的。因為多元高斯分布的任意有限自己還是多遠高斯分布,所以先驗和後驗也都可以都用高斯建模。

實踐中,高斯過程先驗(乙個多元高斯)的均值通常統一設為0(因為我們通常對於需要擬合的函式取值並沒有太多先驗知識),而協方差矩陣則用核函式k(x,x')來建模(i.e., cov[yi, yj] = k(xi, xj) = η exp(-ρ * sed(xi, xj)) )。協方差描述的是,當x變化時,y是如何變化的.而核函式(高斯核)的特點是x,x'越近,返回值越大。用高斯核作協方差可以理解為乙個小的x擾動會導致較小的y變化,乙個大的x變化會導致較大的y變化。

此外,現實中觀測值一定會受到雜訊的影響,因此我們對觀測值建模需要引入一項高斯白雜訊:yi ~ f(xi) +  εi. ε ~ n(0, σ)。這樣協方差函式就改寫成 cov[yi, yj] = k(xi, xj) + σ δij.

注意高斯過程的先驗建模中有引數θ = , 但為了強調高斯過程是乙個非參方法,這些引數稱為「超參」。這些引數將用貝葉斯推斷來估計。

總結一下,我們有各種假設先驗(θ的先驗分布,f(x) ~ n(0, k(x, x'))),高斯似然 y | θ ~ n(0, k(x, x) + σ i), 然後根據y的觀測值可以推斷出超參的分布。**就更簡單了,高斯過程的乙個優點就是後驗分布可以直接求得解析解,要**的點記為(x*, f(x*)), 後驗分布可得 f(x*) | x, x*, y ~ n(μ*, σ*), 其中 μ* = k(x*, x) (k(x, x) + σ i)^(-1) y, σ* = k(x*, x*) - k(x*, x) (k(x, x) + σ i)^(-1) k(x, x*).詳細推導可見gaussian processes for machine learning一書的sec. 2.2, 此外可以記住的乙個常見結論如下:

然後是**實踐,用20個帶噪的正弦函式觀測點來擬合函式:

#%matplotlib inline

import pymc3 as pm

import numpy as np

import theano.tensor as tt

import matplotlib.pyplot as plt

from scipy.spatial.distance import cdist

if __name__ == "__main__":

np.random.seed(1)

squared_distance = lambda x, y: cdist(x.reshape(-1,1), y.reshape(-1,1)) ** 2 #sed function

n = 20 # number of training points.

n = 100 # number of test points.

np.random.seed(1)

f = lambda x: np.sin(x).flatten()

x = np.random.uniform(0, 10, size=n)

y = np.random.normal(np.sin(x), np.sqrt(0.01))

plt.plot(x, y, 'o')

plt.xlabel('$x$', fontsize=16)

plt.ylabel('$f(x)$', fontsize=16, rotation=0)

with pm.model() as gp:

mu = np.zeros(n)

eta = pm.halfcauchy('eta', 0.1)

rho = pm.halfcauchy('rho', 1)

sigma = pm.halfcauchy('sigma', 1)

d = squared_distance(x, x) #sed(x,x)

k = tt.fill_diagonal(eta * pm.math.exp(-rho * d), eta + sigma) #(k(x, x) + σ i)

obs = pm.mvnormal('obs', mu, cov=k, observed=y)

test_points = np.linspace(0, 10, 100)

d_pred = squared_distance(test_points, test_points) #sed(x*,x*)

d_off_diag = squared_distance(x, test_points) #sed(x,x*) n * n

k_oo = eta * pm.math.exp(-rho * d_pred) #k(x*,x*)

k_o = eta * pm.math.exp(-rho * d_off_diag) #k(x,x*)

inv_k = tt.nlinalg.matrix_inverse(k)

mu_post = pm.deterministic('mu_post', pm.math.dot(pm.math.dot(k_o.t, inv_k), y))

sigma_post = pm.deterministic('sigma_post', k_oo - pm.math.dot(pm.math.dot(k_o.t, inv_k), k_o))

step = pm.metropolis()

start = pm.find_map()

trace = pm.sample(1000, step = step, start=start, nchains = 1)

varnames = ['eta', 'rho', 'sigma']

chain = trace[100:]

pm.traceplot(chain, varnames)

plt.figure()

y_pred = [np.random.multivariate_normal(m, s) for m,s in zip(chain['mu_post'], chain['sigma_post'])]

for yp in y_pred:

plt.plot(test_points, yp, 'r-', alpha=0.1)

plt.plot(x, y, 'bo')

plt.xlabel('$x$', fontsize=16)

plt.ylabel('$f(x)$', fontsize=16, rotation=0)

輸出:

高斯樸素貝葉斯

在前面的幾篇部落格中,對樸素貝葉斯的理論知識進行了乙個學習與總結,接下來希望對sklearn庫中的樸素貝葉斯分類器作進一步的學習和說明。bayes.gaussiannb priors none,var smoothing 1e 09 包含兩個引數 prior 表示類的先驗概率 即,沒有條件下的p y...

貝葉斯分析

1先來說一下貝葉斯統計與經典統計的不同之處 簡單說,頻率派認為估計物件 引數 是乙個未知的固定值。而貝葉斯卻認為未知的引數都是隨機變數。我曾經見到這麼個不錯的例子 我們要通過一些事實估計 愛因斯坦在1905年12月25日晚上八點吸菸 的真假。定義引數 那麼頻率派認為,愛因斯坦有沒有曾經在這時刻吸菸是...

貝葉斯 01 初識貝葉斯

分割線 分割線 分割線 分割線 分割線 分割線 分割線 分割線 分割線 分割線 分割線 分割線 分割線 最先知道貝葉斯公式還是四年前的概率論和數理統計課上,時間也很久了,具體內容早已經忘記,不過畢竟曾經學過,重新看過還是得心應手的。大概用兩三篇的內容來介紹一下貝葉斯,以及機器學習中很重要的一部分 樸...