Gaussian Process Modeling - Revision history

Admin at 17:32, 24 June 2026

2026-06-24T17:32:23Z

← Older revision		Revision as of 17:32, 24 June 2026
Line 1:		Line 1:
	GPM = Gaussian Process Modeling~~. It's the~~ rare statistical tool that hands you uncertainty estimates as a built-in feature - no extra computation, no cross-validation headaches, the math just coughs them up as part of the deal. This page is a guided tour through doing Gaussian Process Modeling in Python, from hand-rolled covariance functions all the way to full Bayesian MCMC sampling.		GPM = Gaussian Process Modeling

			The rare statistical tool that hands you uncertainty estimates as a built-in feature - no extra computation, no cross-validation headaches, the math just coughs them up as part of the deal.

			This page is a guided tour through doing Gaussian Process Modeling in Python, from hand-rolled covariance functions all the way to full Bayesian MCMC sampling.

	== Why GPs Deserve the Hype ==		== Why GPs Deserve the Hype ==

Admin: Admin moved page GPM to Gaussian Process Modeling

2026-06-24T17:31:40Z

Admin moved page GPM to Gaussian Process Modeling

← Older revision

Revision as of 17:31, 24 June 2026

Admin: Replace all em-dashes (—) with regular dashes (-) (via update-page on MediaWiki MCP Server)

2026-06-24T17:31:39Z

Replace all em-dashes (—) with regular dashes (-) (via update-page on MediaWiki MCP Server)

Show changes

Admin: /* Where This Page Lives */

2026-06-24T17:30:27Z

Where This Page Lives

← Older revision		Revision as of 17:30, 24 June 2026
Line 370:		Line 370:

	The field is actively working on making these approximations better and cheaper, so expect this landscape to keep shifting.		The field is actively working on making these approximations better and cheaper, so expect this landscape to keep shifting.

	~~== Where This Page Lives ==~~

	[[Category:Data Analysis]]		[[Category:Data Analysis]]
	[[Category:Python]]		[[Category:Python]]
	[[Category:Monte Carlo]]		[[Category:Monte Carlo]]

Admin: gave this page the Ballmer treatment — better headings, actual personality, fixed bare URLs and typos (via update-page on MediaWiki MCP Server)

2026-06-24T17:28:04Z

gave this page the Ballmer treatment — better headings, actual personality, fixed bare URLs and typos (via update-page on MediaWiki MCP Server)

Show changes

Admin: /* Flags */

2017-11-06T11:50:47Z

Flags

← Older revision		Revision as of 11:50, 6 November 2017
Line 307:		Line 307:

	This is useful when the likelihood function or model is unusual or difficult to fit using the gradient ascent optimization methods built into scikit-learn.		This is useful when the likelihood function or model is unusual or difficult to fit using the gradient ascent optimization methods built into scikit-learn.

			===PyMC3===

			[[PyMC]]3 is a Python module for probabilistic programming for fitting a Bayesian model to data. PyMC3 uses a Theano backend (analogous to GPflow using TensorFlow as the backend). It provides functions and objects for specifying covariance and prior distribution kernels. Like GPflow, we must specify these as tensor variables.

			To import it, we import PyMC3 and Theano:

			<pre>
			import pymc3 as pm
			import theano.tensor as tt
			</pre>

			A Model() context is needed to declare parameters, variables, and functions within a Bayesian model. (Se PyMC3 documentation [http://docs.pymc.io/ here]). For example, we'll construct a Matern distribution again:

			<pre>
			with pm.Model() as gp_fit:
			rho = pm.Gamma('rho', 1, 1)
			eta = pm.Gamma('eta', 1, 1)
			k = eta * pm.gp.cov.Matern32(1, rho)
			</pre>

			Next, the mean function is specified (this will be zero if none is specified), and a half-Cauchy distribution is specified for a noise variable:

			<pre>
			with gp_fit:
			M = pm.gp.mean.Zero()
			sigma = pm.HalfCauchy('sigma', 2.5)
			</pre>

			PyMC implements Gaussian process models using the GP class. This requires a mean function, covariance function, and observation error function. These have been specified above.

			The outcomes of the GP have been observed: the data are contained in the X and y variables. We provide these to the GP object that we create:

			<pre>
			with gp_fit:
			y_obs = pm.gp.GP('y_obs', mean_func=M, cov_func=K, sigma=sigma, observed={'X':X, 'Y':y})
			</pre>

			To fit the model, we call the sample() function from within the Model context. By default, PyMC uses the No U-turn Sampler (NUTS), which is a version of Hamiltonian Monte Carlo (which, in turn, is a form of Markov chain Monte Carlo). Like GPflow, this uses Theano's internals to compute gradients for arbitrary models.

			<pre>
			with gp_fit:
			trace = pm.sample(2000, n_init = 20000)

			pm.traceplot( trace[-1000:], varnames=['rho', 'sigma', 'eta'])
			</pre>

			PyMC allows predictive sampling using the fit model. The sample_gp function implements sampling over a grid of points:

			<pre>
			Z = np.linspace(-6, 6, 100).reshape(-1, 1)

			with gp_fit:
			gp_samples = pm.gp.sample_gp( trace[1000:], y_obs, Z, samples=50)
			</pre>

			This will generate 50 different random, Monte Carlo models that fit the data points.

			Why use PyMC? MCMC is costly for very large data sets (because the log-probability must be evaluated at each iteration). Variational inference methods approximate the posterior with a cheaper and simpler approximation. While this makes it less accurate, it is appropriate for some problems.

			Future work/methods will focus on improving the quality and lowering the cost of the approximation.

			Variational Gaussian Approximation is implemented in GPflow.

			Automatic Differentiation Variational Inference is implemented in PyMC.

	=Flags=		=Flags=

Admin: /* Python Libraries */

2017-11-06T10:50:48Z

Python Libraries

← Older revision		Revision as of 10:50, 6 November 2017
Line 13:		Line 13:

	These packages principally implement covariance functions that can be used to describe non-linearity in data, and methods for fitting GPM parameters.		These packages principally implement covariance functions that can be used to describe non-linearity in data, and methods for fitting GPM parameters.

			===Sample Data===

			To create sample data:

			<pre>
			import numpy as np

			def exponential_cov(x, y, params):
			return params[0] * np.exp( -0.5 * params[1] * np.subtract.outer(x, y)**2)

			def conditional(x_new, x, y, params):
			B = exponential_cov(x_new, x, params)
			C = exponential_cov(x, x, params)
			A = exponential_cov(x_new, x_new, params)

			mu = np.linalg.inv(C).dot(B.T).T.dot(y)
			sigma = A - B.dot(np.linalg.inv(C).dot(B.T))

			return(mu.squeeze(), sigma.squeeze())
			</pre>

			We'll create a Gaussian process prior with hyperparameters theta0 = 1 and theta1 = 10, along with a zero mean:

			<pre>
			import matplotlib.pylab as plt

			theta = [1, 10]
			sigma_0 = exponential_cov(0, 0, theta)
			xpts = np.arange(-3, 3, step=0.01)
			plt.errorbar(xpts, np.zeros(len(xpts)), yerr=sigma_0, capsize=0)
			</pre>

			Now, pick an arbitrary starting point:

			<pre>
			x = [1.]
			y = [np.random.normal(scale=sigma_0)]
			</pre>

			Now update the confidence band to account for the new sample point by using the covariance function to generate point-wise intervals:

			<pre>
			sigma_1 = exponential_cov(x, x, theta)

			def predict(x, data, kernel, params, sigma, t):
			k = [kernel(x, y, params) for y in data]
			Sinv = np.linalg.inv(sigma)
			y_pred = np.dot(k, Sinv).dot(t)
			sigma_new = kernel(x, x, params) - np.dot(k, Sinv).dot(k)
			return y_pred, sigma_new

			x_pred = np.linspace(-3, 3, 1000)
			predictions = [predict(i, x, exponential_cov, theta, sigma_1, y) for i in x_pred]
			</pre>

			To plot:

			<pre>
			y_pred, sigmas = np.transpose(predictions)
			plt.errorbar(x_pred, y_pred, yerr=sigmas, capsize=0)
			plt.plot(x, y, "ro")
			</pre>

			Another point can be added as well:

			<pre>
			m, s = conditional([-0.7], x, y, θ)
			y2 = np.random.normal(m, s)
			</pre>

			To update the covariance to account for the new point:

			<pre>
			x.append(-0.7)
			y.append(y2)

			sigma_2 = exponential_cov(x, x, theta)
			predictions = [predict(i, x, exponential_cov, theta, sigma_2, y) for i in x_pred]
			</pre>

			To sample multiple points at once:

			<pre>
			x_more = [-2.1, -1.5, 0.3, 1.8, 2.5]
			mu, s = conditional(x_more, x, y, theta)
			y_more = np.random.multivariate_normal(mu, s)
			</pre>

			Now the covariance can be updated again:

			<pre>
			x += x_more
			y += y_more.tolist()

			σ_new = exponential_cov(x, x, theta)
			predictions = [predict(i, x, exponential_cov, theta, sigma_new, y) for i in x_pred]

			y_pred, sigmas = np.transpose(predictions)
			plt.errorbar(x_pred, y_pred, yerr=sigmas, capsize=0)
			plt.plot(x, y, "ro")
			</pre>

	===Scikit-Learn===		===Scikit-Learn===

Admin: /* GPFlow */

2017-11-06T10:43:27Z

GPFlow

@@ Line 120: / Line 120: @@
 <pre>
 k = GPflow.kernels.Matern32(1, variance=1, lengthscales=1.2)
-m = GPflow.gpr.GPR(X, Y, kern=k)
+Now, to fit the model, call optimize():
 print(m)
 </pre>
 =Flags=

Admin: /* GPFlow */

2017-11-06T10:09:28Z

GPFlow

@@ Line 93: / Line 93: @@
 sigma_0 = exponential_cov(0, 0, theta)
 y = [np.random.normal(scale=sigma_0)]
 Y = y.reshape(-1,1)
 </pre>
@@ Line 102: / Line 104: @@
 def exponential_cov(x, y, params):
      return params[0] * np.exp(-0.5*params[1]*np.subtract.outer(x, y)**2)
 </pre>

Admin: /* GPFlow */

2017-11-06T09:51:08Z

GPFlow

← Older revision		Revision as of 09:51, 6 November 2017
Line 101:		Line 101:
	import numpy as np		import numpy as np
	def exponential_cov(x, y, params):		def exponential_cov(x, y, params):
	return params[0] * np.exp( -0.5 * params[1] * np.subtract.outer(x, y)**2)		return params[0] * np.exp(-0.5params[1]np.subtract.outer(x, y)**2)
	</pre>		</pre>