# Generate the randomized sampleX = 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 = 30w0 = 2.0
w1 = 4.5x = 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 = 30x = 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 distributionmus = [[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 = -3N = 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 distributionmus = [0, 0]
covs = [[1, 0.8], [0.8, 1]]fig, ax = plt.subplots(1, 1)max_val = 3
min_val = -3N = 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 npdef plot_gpr_samples(gpr_model, n_samples, ax):"""Borrow from https://scikit-learn.org/stable/auto_examples/gaussian_process/plot_gpr_prior_posterior.htmlPlot samples drawn from the Gaussian process model.If the Gaussian process model is not trained then the drawn samples aredrawn from the prior distribution. Otherwise, the samples are drawn fromthe posterior distribution. Be aware that a sample here corresponds to afunction.Parameters----------gpr_model : `GaussianProcessRegressor`A :class:`~sklearn.gaussian_process.GaussianProcessRegressor` model.n_samples : intThe number of samples to draw from the Gaussian process distribution.ax : matplotlib axisThe 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 kernelkernel = (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 logreturn 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 kernelkernel = 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 sampleX = 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 visualizationplt.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.int32fig = m.plot()
GPy.plotting.show(fig, filename='basic_gp_regression_notebook')

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[15], line 41 #np.int = np.int323 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 by145     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.figureAttributeError: '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 [ ]: