GPM - Revision history

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>

Admin: /* GPFlow */

2017-11-06T09:36:55Z

GPFlow

← Older revision		Revision as of 09:36, 6 November 2017
Line 82:		Line 82:

	===GPFlow===		===GPFlow===

			GPFlow is software derived from GPy software, a Gaussian Process fitting software from teh Sheffield machine learning group. Both provide a set of classes for creating Gaussian process models, and a library of kernels for mean and covariate functions. Under the hood, GPFlow uses the TensorFlow library. The use of TensorFlow allows fitting larger and more complex Gaussian Process models, and more complicated methods like variational inference.

			The API requires a tabulated input for predictors (features) and outputs.

			To reshape y to a tabular format:

			<pre>
			theta = [1, 10]
			sigma_0 = exponential_cov(0, 0, theta)
			y = [np.random.normal(scale=sigma_0)]
			Y = y.reshape(-1,1)
			</pre>

			(a reminder, the exponential_cov function is defined as:)

			<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)
			</pre>

	=Flags=		=Flags=

Admin: /* Scikit-Learn */

2017-11-06T09:24:08Z

Scikit-Learn

← Older revision		Revision as of 09:24, 6 November 2017
Line 80:		Line 80:
	y_pred, sigma = gp.predict(x_pred, return_std=True)		y_pred, sigma = gp.predict(x_pred, return_std=True)
	</pre>		</pre>

			===GPFlow===

	=Flags=		=Flags=

Admin: /* Scikit-Learn */

2017-11-06T09:19:33Z

Scikit-Learn

@@ Line 50: / Line 50: @@
 </pre>
-Or, to do this all in a single call:
+To do this all in a single call, while also setting the optimizer:
 <pre>
@@ Line 57: / Line 57: @@
      n_restarts_optimizer=0, normalize_y=False,
      optimizer='fmin_l_bfgs_b', random_state=None)
 </pre>

Admin at 08:10, 6 November 2017

2017-11-06T08:10:02Z

← Older revision		Revision as of 08:10, 6 November 2017
Line 58:		Line 58:
	optimizer='fmin_l_bfgs_b', random_state=None)		optimizer='fmin_l_bfgs_b', random_state=None)
	</pre>		</pre>

			=Flags=

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

Admin: /* Scikit-Learn */

2017-11-06T00:18:56Z

Scikit-Learn

← Older revision		Revision as of 00:18, 6 November 2017
Line 38:		Line 38:

	<pre>		<pre>
	kernel = ConstantKernel() + Matern(length_scale=2, nu=3/2) + WhiteKernel(noise_level=1)		kernel = ConstantKernel()
			+ Matern(length_scale=2, nu=3/2)
			+ WhiteKernel(noise_level=1)
	</pre>		</pre>