高斯过程的数学理解

# Generate the randomized sample

X = np.linspace(start=0, stop=10, num=100).reshape(-1, 1)
y = np.squeeze(X * np.sin(X)) + np.random.randn(X.shape[0])
y = y.reshape(-1, 1)

# RBF kernel
kernel = GPy.kern.RBF(input_dim=1, variance=1., lengthscale=1.)

# Gaussian process regression using RBF kernel
m = GPy.models.GPRegression(X, y, kernel)

# Load the diabetes dataset and create a dataframe
diabetes = datasets.load_diabetes()
df = pd.DataFrame(diabetes.data, columns=diabetes.feature_names)

# dataset
X = df[['age', 'sex', 'bmi', 'bp']].values
y = df[['target']].values

# set kernel
kernel = GPy.kern.RBF(input_dim=4, variance=1.0, lengthscale=1.0, ARD = True)

# create GPR model
m = GPy.models.GPRegression(X, y, kernel = kernel) 

Mathematical understanding of Gaussian Process | by Yuki Shizuya | The Quantastic Journal | Jun, 2024 | Medium

七、代码参考

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from scipy.stats import norm, multivariate_normal
from sklearn.linear_model import LinearRegression
from sklearn import datasets

import GPy

Linear data

In [15]:

n_samples = 30

w0 = 2.0
w1 = 4.5

x = np.sort(np.random.rand(n_samples))
y = w1 * x + w0 + np.random.rand(n_samples)

x = x.reshape(-1, 1)
y = y.reshape(-1, 1)

lr = LinearRegression()
lr.fit(x, y)

print('coefficient:', lr.coef_)
print('intercept:', lr.intercept_)

fig, ax = plt.subplots(1, 2, figsize=(10, 5))

ax[0].scatter(x, y, color='blue', label='input data')
ax[0].set_xlabel('x')
ax[0].set_ylabel('y')

ax[1].scatter(x, y, color='blue', label='input data')
ax[1].plot(x, lr.predict(x), color='black', label='linear regression line')
ax[1].set_xlabel('x')
ax[1].set_ylabel('y')
ax[1].legend()
plt.show()

coefficient: [[4.60793678]]
intercept: [2.38216139]

Non-linear data

In [18]:

n_samples = 30

x = np.sort(np.random.rand(n_samples))
y = np.sin(2.0 * np.pi * x) + np.random.rand(n_samples)

x = x.reshape(-1, 1)
y = y.reshape(-1, 1)

lr = LinearRegression()
lr.fit(x, y)

print('coefficient:', lr.coef_)
print('intercept:', lr.intercept_)

fig, ax = plt.subplots(1, 2, figsize=(10, 5))

ax[0].scatter(x, y, color='blue', label='input data')
ax[0].set_xlabel('x')
ax[0].set_ylabel('y')

ax[1].scatter(x, y, color='blue', label='input data')
ax[1].plot(x, lr.predict(x), color='black', label='linear regression line')
ax[1].set_xlabel('x')
ax[1].set_ylabel('y')
ax[1].legend()

plt.show()

coefficient: [[-2.12667242]]
intercept: [1.5927719]

Radial basis function regression

In [22]:

# the mu values for radial basis function
mus = np.array([i for i in range(-4, 4, 1)])

# 1-D input data
x = np.linspace(-5, 5, 500)

# for visualization
fig, ax = plt.subplots(1, 1)

# calculate the pdf of each mu
for mu in mus:
    rv = norm(loc=mu, scale=1.0)
    ax.plot(x, rv.pdf(x), color='blue')
    
ax.set_xlabel('x')
ax.set_title('Radial basis function (1 dimension)')
plt.show()

In [33]:

# the mu values for radial basis function
mus = np.array([i for i in range(-4, 5, 1)])

# weights for the spline function
weights = np.array([-0.48, -0.52, 0.1, -0.9, 0.23, 0.57, 0.28, -0.69, 0.46])

# 1-D input data
x = np.linspace(-5, 5, 500)

y = np.zeros(x.shape)

for i in range(len(x)):
    for j in range(len(weights)):
        y[i] += weights[j] * np.exp(-(x[i]-mus[j])**2)

# for visualization
fig, ax = plt.subplots(1, 1)
ax.plot(x, y, color='black')
ax.set_xlabel('x')
ax.set_title('Radial basis function regression')
plt.grid()
plt.show()

Merginal distribution visualization

In [25]:

# Marginal distribution

mus = [[0, 0] for _ in range(3)]
covs = [[[1, 0.2], [0.2, 1]], [[1, 0.5], [0.5, 1]], [[1, 0.8], [0.8, 1]]]

fig, ax = plt.subplots(1, 3, figsize=(15, 5))

max_val = 3
min_val = -3

N = 200
x = np.linspace(min_val, max_val, N)
y = np.linspace(min_val, max_val, N)
X, Y = np.meshgrid(x, y)

pos  = np.dstack((X, Y))

for i in range(3):
    mu = mus[i]
    cov = covs[i]
    
    rv   = multivariate_normal(mu, cov)
    Z    = rv.pdf(pos)

    ax[i].contour(X, Y, Z)
    ax[i].set_aspect('equal', adjustable='box')
    
    rv1 = norm(mu[0], cov[0][0])
    z1 = rv1.pdf(x)
    ax[i].plot(x, z1-max_val)
    
    rv2 = norm(mu[1], cov[1][1])
    z2 = rv2.pdf(y)
    ax[i].plot(z2-max_val, y)
    
    ax[i].set_xticks([])
    ax[i].set_yticks([])

plt.grid()
plt.show()

Conditional distribution visualization

In [40]:

# Conditional distribution

mus = [0, 0]
covs = [[1, 0.8], [0.8, 1]]

fig, ax = plt.subplots(1, 1)

max_val = 3
min_val = -3

N = 200
x = np.linspace(min_val, max_val, N)
y = np.linspace(min_val, max_val, N)
X, Y = np.meshgrid(x, y)

pos  = np.dstack((X, Y))
rv   = multivariate_normal(mu, cov)
Z    = rv.pdf(pos)

ax.contour(X, Y, Z)
ax.set_aspect('equal', adjustable='box')

for c in [-1.5, 0, 1.5]:
    cond_mu = mu[1] + covs[1][0] * (1 / covs[0][0]) * (c - mu[0])
    cond_cov = covs[1][1] - covs[1][0] * (1 / covs[0][0]) * covs[0][1]

    rv1 = norm(cond_mu, cond_cov)
    z1 = rv1.pdf(x)
    ax.plot(z1+c, x, '--', color='black')
    ax.plot([c for _ in range(len(y))], y, color='black')

plt.grid()
plt.show()

Sampling from Gaussain process using various kernel funcions

In [42]:

import matplotlib.pyplot as plt
import numpy as np


def plot_gpr_samples(gpr_model, n_samples, ax):
    """
    Borrow from https://scikit-learn.org/stable/auto_examples/gaussian_process/plot_gpr_prior_posterior.html
    
    Plot samples drawn from the Gaussian process model.

    If the Gaussian process model is not trained then the drawn samples are
    drawn from the prior distribution. Otherwise, the samples are drawn from
    the posterior distribution. Be aware that a sample here corresponds to a
    function.

    Parameters
    ----------
    gpr_model : `GaussianProcessRegressor`
        A :class:`~sklearn.gaussian_process.GaussianProcessRegressor` model.
    n_samples : int
        The number of samples to draw from the Gaussian process distribution.
    ax : matplotlib axis
        The matplotlib axis where to plot the samples.
    """
    
    x = np.linspace(0, 5, 100)
    X = x.reshape(-1, 1)

    y_mean, y_std = gpr_model.predict(X, return_std=True)
    y_samples = gpr_model.sample_y(X, n_samples)

    for idx, single_prior in enumerate(y_samples.T):
        ax.plot(
            x,
            single_prior,
            linestyle="--",
            alpha=0.7,
            label=f"Sampled function #{idx + 1}",
        )
    ax.plot(x, y_mean, color="black", label="Mean")
    ax.fill_between(
        x,
        y_mean - y_std,
        y_mean + y_std,
        alpha=0.1,
        color="black",
        label=r"$\pm$ 1 std. dev.",
    )
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_ylim([-3, 3])

Example dataset to visualize sampling values from various kernels

In [56]:

rng = np.random.RandomState(4)
X_train = rng.uniform(0, 5, 10).reshape(-1, 1)
y_train = np.sin((X_train[:, 0] - 2.5) ** 2)
n_samples = 5

In [63]:

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel, DotProduct, ExpSineSquared

In [59]:

# Linear kernel

kernel = (
    ConstantKernel(0.1, (0.01, 10.0)) + DotProduct(sigma_0=0)
)
gpr = GaussianProcessRegressor(kernel=kernel, random_state=0)

fig, axs = plt.subplots(nrows=2, sharex=True, sharey=True, figsize=(10, 8))

# plot prior
plot_gpr_samples(gpr, n_samples=n_samples, ax=axs[0])
axs[0].set_title("Samples from Gaussian process using linear kernel")

# plot posterior
gpr.fit(X_train, y_train)
plot_gpr_samples(gpr, n_samples=n_samples, ax=axs[1])
axs[1].scatter(X_train[:, 0], y_train, color="red", zorder=10, label="Observations")
axs[1].legend(bbox_to_anchor=(1.05, 1.5), loc="upper left")
axs[1].set_title("Samples from posterior distribution")

fig.suptitle("Dot-product kernel", fontsize=18)
plt.tight_layout()

/opt/homebrew/anaconda3/envs/smoothing-env/lib/python3.8/site-packages/sklearn/gaussian_process/kernels.py:284: RuntimeWarning: divide by zero encountered in log
  return np.log(np.hstack(theta))
/opt/homebrew/anaconda3/envs/smoothing-env/lib/python3.8/site-packages/sklearn/gaussian_process/_gpr.py:659: ConvergenceWarning: lbfgs failed to converge (status=2):
ABNORMAL_TERMINATION_IN_LNSRCH.

Increase the number of iterations (max_iter) or scale the data as shown in:
    https://scikit-learn.org/stable/modules/preprocessing.html
  _check_optimize_result("lbfgs", opt_res)
/opt/homebrew/anaconda3/envs/smoothing-env/lib/python3.8/site-packages/sklearn/gaussian_process/kernels.py:419: ConvergenceWarning: The optimal value found for dimension 0 of parameter k2__sigma_0 is close to the specified lower bound 1e-05. Decreasing the bound and calling fit again may find a better value.
  warnings.warn(

In [61]:

# Gaussian kernel

kernel = 1.0 * RBF(length_scale=1.0, length_scale_bounds=(1e-1, 10.0))
gpr = GaussianProcessRegressor(kernel=kernel, random_state=0)

fig, axs = plt.subplots(nrows=2, sharex=True, sharey=True, figsize=(10, 8))

# plot prior
plot_gpr_samples(gpr, n_samples=n_samples, ax=axs[0])
axs[0].set_title("Samples from Gaussian process using Gaussian kernel")

# plot posterior
gpr.fit(X_train, y_train)
plot_gpr_samples(gpr, n_samples=n_samples, ax=axs[1])
axs[1].scatter(X_train[:, 0], y_train, color="red", zorder=10, label="Observations")
axs[1].legend(bbox_to_anchor=(1.05, 1.5), loc="upper left")
axs[1].set_title("Samples from posterior distribution")

fig.suptitle("Radial Basis Function kernel", fontsize=18)
plt.tight_layout()

In [65]:

kernel = 1.0 * ExpSineSquared(
    length_scale=1.0,
    periodicity=1.0,
    length_scale_bounds=(0.1, 10.0),
    periodicity_bounds=(1.0, 10.0),
)
gpr = GaussianProcessRegressor(kernel=kernel, random_state=0)

fig, axs = plt.subplots(nrows=2, sharex=True, sharey=True, figsize=(10, 8))

# plot prior
plot_gpr_samples(gpr, n_samples=n_samples, ax=axs[0])
axs[0].set_title("Samples from Gaussian process using periodic kernel")

# plot posterior
gpr.fit(X_train, y_train)
plot_gpr_samples(gpr, n_samples=n_samples, ax=axs[1])
axs[1].scatter(X_train[:, 0], y_train, color="red", zorder=10, label="Observations")
axs[1].legend(bbox_to_anchor=(1.05, 1.5), loc="upper left")
axs[1].set_title("Samples from posterior distribution")

fig.suptitle("Exp-Sine-Squared kernel", fontsize=18)
plt.tight_layout()

/opt/homebrew/anaconda3/envs/smoothing-env/lib/python3.8/site-packages/sklearn/gaussian_process/kernels.py:429: ConvergenceWarning: The optimal value found for dimension 0 of parameter k2__periodicity is close to the specified upper bound 10.0. Increasing the bound and calling fit again may find a better value.
  warnings.warn(

Example: Gaissian process model for one-dimensional data

The code heavily borrowed from Jupyter Notebook Viewer

In [11]:

# Generate the randomized sample

X = np.linspace(start=0, stop=10, num=100).reshape(-1, 1)
y = np.squeeze(X * np.sin(X)) + np.random.randn(X.shape[0])
y = y.reshape(-1, 1)

In [12]:

# The input data's visualization

plt.scatter(X, y)
plt.show()

In [3]:

# RBF kernel
kernel = GPy.kern.RBF(input_dim=1, variance=1., lengthscale=1.)

In [13]:

# Gaussian process regression using RBF kernel
m = GPy.models.GPRegression(X, y, kernel)

In [14]:

from IPython.display import display
display(m)

Model: GP regression
Objective: 204.66508168561933
Number of Parameters: 3
Number of Optimization Parameters: 3
Updates: True

In [15]:

#np.int = np.int32

fig = m.plot()
GPy.plotting.show(fig, filename='basic_gp_regression_notebook')

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[15], line 4
      1 #np.int = np.int32
      3 fig = m.plot()
----> 4 GPy.plotting.show(fig, filename='basic_gp_regression_notebook')

File /opt/homebrew/anaconda3/envs/smoothing-env/lib/python3.8/site-packages/GPy/plotting/__init__.py:150, in show(figure, **kwargs)
    142 def show(figure, **kwargs):
    143     """
    144     Show the specific plotting library figure, returned by
    145     add_to_canvas().
   (...)
    148     for showing/drawing a figure.
    149     """
--> 150     return plotting_library().show_canvas(figure, **kwargs)

File /opt/homebrew/anaconda3/envs/smoothing-env/lib/python3.8/site-packages/GPy/plotting/matplot_dep/plot_definitions.py:96, in MatplotlibPlots.show_canvas(self, ax, **kwargs)
     95 def show_canvas(self, ax, **kwargs):
---> 96     ax.figure.canvas.draw()
     97     return ax.figure

AttributeError: 'dict' object has no attribute 'figure'

In [ ]:

Example : Gaussian process model for multiple dimension data

In [21]:

# Load the diabetes dataset and create a dataframe
diabetes = datasets.load_diabetes()
df = pd.DataFrame(diabetes.data, columns=diabetes.feature_names)
df.head()

Out[21]:

In [23]:

# Add the target variable to the dataframe
df['target'] = diabetes.target

In [36]:

# training dataset
X = df[['age', 'sex', 'bmi', 'bp']].values
y = df[['target']].values

In [37]:

# set kernel
kernel = GPy.kern.RBF(input_dim=4, variance=1.0, lengthscale=1.0, ARD = True)

In [38]:

# create GPR model
m = GPy.models.GPRegression(X, y, kernel = kernel, normalizer = True) 

In [39]:

display(m)

Model: GP regression
Objective: 580.753895481785
Number of Parameters: 6
Number of Optimization Parameters: 6
Updates: True

In [46]:

intervals = m.predict_quantiles(X)

In [47]:

df['quantile_2.5%'] = intervals[0]
df['quantile_97.5%'] = intervals[1]

In [48]:

compare_dataset = df[['age', 'sex', 'bmi', 'bp', 'target', 'quantile_2.5%', 'quantile_97.5%']]

In [49]:

compare_dataset.head(20)

Out[49]:

In [ ]: