"""New Module to implement tasks relating to Gaussian Processes.
2013-03-14 18:40 IJMC: Begun.
"""
import numpy as np
[docs]def squaredExponentialKernel(params, x, x1=None):
"""Construct a squared-exponential kernel with specified perameters.
:INPUTS:
x : 1D NumPy array
Input independent data.
params : sequence
params[0] -- h : scalar
Amplitude scaling
params[1] -- L : scalar
Length scaling
params[2] -- s : OPTIONAL scalar
White noise component from :func:`whiteNoiseKernel` (optional!)
x1 : 1D NumPy Array
Dependent data (e.g.)
:NOTES:
Computes k(x, x') = h^2 exp(- [ (x - x') / 2L ]^2 )
:RETURNS:
k, the kernel matrix.
:REFERENCE:
Roberts et al. 2012, Eq. 13.
:EXAMPLE:
::
import gaussianprocess as gp
import numpy as np
import pylab as py
x0 = np.arange(50.) # independent data
k = gp.squaredExponentialKernel([1, 5], x0)
sample_draw = np.random.multivariate_normal(np.zeros(x0.size), k, 1).ravel()
"""
# 2013-03-14 18:43 IJMC: Created.
if x1 is None:
x1 = x
h, L = params[0:2]
if h<0:
h = np.abs(h)
if L<0:
L = np.abs(L)
ret = h**2 * np.exp(-0.5 * ((x.reshape(x.size, 1) - x1.reshape(1, x1.size)) / L)**2)
if len(params)>=3:
ret += whiteNoiseKernel(params[2], x.size)
return ret
[docs]def whiteNoiseKernel(s, x):
"""Construct a white-noise kernel (diagonal matrix) with specified sigma.
:INPUTS:
x : 1D NumPy array, or scalar
Input independent data, or size of desired matrix.
s : scalar
Sigma [i.e., sqrt(variance) ] of the desired white noise.
:NOTES:
Computes k(i, j) = delta_{ij} s^2
:RETURNS:
k, the kernel matrix.
:REFERENCE:
Roberts et al. 2012, Eq. 12.
:EXAMPLE:
::
import gaussianprocess as gp
import numpy as np
import pylab as py
x0 = np.arange(50.) # independent data
k = gp.whiteNoiseKernel(1, x0)
sample_draw = np.random.multivariate_normal(np.zeros(x.size), k, 1).ravel()
"""
# 2013-03-14 18:43 IJMC: Created.
if s < 0:
s = np.abs(s)
if hasattr(x, '__iter__') and x.size>0:
ret = np.diag(np.ones(x.size) * s**2)
else:
ret = np.diag(np.ones(x) * s**2)
return ret
[docs]def negLogLikelihood(*arg, **kw):
"""Returns negative log likelihood, for minimization routines.
See :func:`logLikelihood` for syntax and details."""
# 2013-03-15 15:02 IJMC: Created
return -logLikelihood(*arg, **kw)
[docs]def logLikelihood(*arg, **kw):
"""Compute log likelihood using Gaussian Process techniques.
:INPUTS:
(fitparams, function, arg1, arg2, ... , depvar, cov_func, cov_args, ncov_params)
OR:
(fitparams, function, arg1, arg2, ... , depvar, cov_func, cov_args, ncov_params, kw)
OR:
(allparams, (args1, args2, ..), npars=(npar1, npar2, ...))
where allparams is an array concatenation of each functions
input parameters.
fitparams : sequence
Parameters used to compute the likelihood. The first values
(length 'ncov_params') are used to compute the covariance
matrix via 'cov_func', and the subsequent values are passed to
'function'.
cov_func : function
Function to generate covariance matrix
(e.g. :func:`squaredExponentialKernel`).
cov_args : tuple
Arguments to be passed to 'cov_func'.
ncov_params: scalar
Number of N parameters used to call 'cov_func'. These
parameters must be the first N values of 'fitparams'!
Computes the likelihood using full covariance matrices, i.e.:
(2 pi |C|)^-0.5 * exp(-0.5 * r.T * C^-1 * r)
where r is the residual vector
(depvar - function(fitparams, arg1, ...)).ravel()
If
'function' need not output a 1D vector, but its output must be
the same shape as 'depvar'. The difference will be
np.ravel()'ed, and the result of the ravel() operation must
be a vector with the same lengths as the 1-D size of
'covariance'.
:OPTIONS:
gaussprior, uniformprior: must be the same length as "params",
which is the concatenation of the paramaters passed to
'cov_func' and 'function'. Elements must be either 'None' or
2-tuples of (mean, std.dev.).
ngaussprior: must be the same length as "params", which is the
concatenation of the parameters passed to 'cov_func' and
'function'! Elements must be either 'None' or 3-tuples of
form: (indices, mean_vector, covariance_matrix).
jointpars : also valid
:SEE_ALSO:
:func:`phasecurves.devfunc` (for a frequentist version without
all these extra covariance-matrix hyperparameters).
"""
# 2013-03-15 11:14 IJMC: Created, from my old function 'devfunc'
# 2014-07-14 15:00 IJMC: Now use numpy.linalg.slogdet instead of .det
if isinstance(arg[-1], dict):
# Surreptiously setting keyword arguments:
kw_func = arg[-1]
kw.update(kw_func)
arg = arg[0:-1]
else:
kw_func = None
if len(arg)==2:
likelihood = likelihood(*arg, **kw)
else:
params = np.array(arg[0], copy=False)
# Put all prior inputs standard form:
if 'gaussprior' in kw and kw['gaussprior'] is not None:
# If any priors are None, redefine them:
gaussprior = []
for pair in kw['gaussprior']:
if pair is None:
gaussprior.append([0, np.inf])
else:
gaussprior.append(pair)
else:
gaussprior = None
if 'uniformprior' in kw and kw['uniformprior'] is not None:
# If any priors are None, redefine them:
uniformprior = []
for pair in kw['uniformprior']:
if pair is None:
uniformprior.append([-np.inf, np.inf])
else:
uniformprior.append(pair)
else:
uniformprior = None
if 'ngaussprior' in kw and kw['ngaussprior'] is not None:
# If any priors are None, redefine them:
ngaussprior = []
for triplet in kw['ngaussprior']:
if triplet is not None and len(triplet)==3:
ngaussprior.append(triplet)
else:
ngaussprior = None
if 'npars' in kw:# simultaneous multiple-function fitting:
ret = 0.
# Excise "npars" kw for recursive calling:
lower_kw = kw.copy()
npars = lower_kw.pop('npars')
# Keep fixed pairs of joint parameters:
if 'jointpars' in kw:
jointpars = kw['jointpars']
for jointpar in jointpars:
params[jointpar[1]] = params[jointpar[0]]
for ii in xrange(len(npars)):
i0 = sum(npars[0:ii])
i1 = i0 + npars[ii]
these_params = arg[0][i0:i1]
ret += likelihood(these_params, *arg[ii+1], **lower_kw)
return ret
else: # Single function-fitting
# Parse these inputs:
function = arg[1]
helperargs = arg[2:len(arg)-4]
depvar = arg[-4]
cov_func, cov_args, ncov_params = arg[-3:]
if not isinstance(cov_args, tuple):
print "Error: cov_args should be a tuple!"
cov_args = (cov_args,)
# (fitparams, function, arg1, arg2, ... , depvar, cov_func, cov_args, ncov_params, kw)
# Compute the residuals:
N = depvar.size
if function is None:
residuals = depvar.reshape(1, N)
else:
residuals = (depvar - function(*((params[ncov_params:],)+helperargs))).reshape(1, N)
# Compute the covariance matrix:
covariance = cov_func(params[0:ncov_params], *cov_args)
# Compute the likelihood:
#pdb.set_trace()
likelihood = -0.5 * (np.dot(residuals, np.dot(np.linalg.inv(covariance), residuals.ravel())) + np.linalg.slogdet(covariance)[1] + N*np.log(2*np.pi))
#if not np.isfinite(likelihood): pdb.set_trace()
# Compute 1D and N-D gaussian, and uniform, prior penalties:
additionalChisq = 0.
if gaussprior is not None:
for param0, gprior in zip(params, gaussprior):
additionalChisq += ((param0 - gprior[0])/gprior[1])**2
if ngaussprior is not None:
for ind, mu, cov in ngaussprior:
dvec = params[ind] - mu
additionalChisq += \
np.dot(dvec.transpose(), np.dot(np.linalg.inv(cov), dvec))
if uniformprior is not None:
for param0, uprior in zip(params, uniformprior):
if (param0 < uprior[0]) or (param0 > uprior[1]):
additionalChisq += 1e16
# L = exp(-0.5 * chi^2), so:
likelihood += -0.5*additionalChisq
#print '%1.3f, %1.2f' % (likelihood, additionalChisq)
# Scale up the residuals so as to impose priors in chi-squared
# space:
return likelihood[0]