Kernel Building¶
The class of kernels that can be used in celerite are sums of terms of the following form:
At a basic level, the interface involves providing the solver with lists of \(\alpha_j\) and \(\beta_j\). It’s useful to know if \(b_j = 0\) and \(d_j = 0\) for any terms because we can use that fact to gain computational efficiency. Therefore, there are 6 arrays that should be provided to the code:
a_real
: an array of \(a_j\) for the terms where \(b_j = 0\) and \(d_j = 0\),c_real
: an array of \(c_j\) for the terms where \(b_j = 0\) and \(d_j = 0\),a_comp
: an array of \(a_j\) for the other terms,b_comp
: an array of \(b_j\) for the other terms,c_comp
: an array of \(c_j\) for the other terms, andd_comp
: an array of \(d_j\) for the other terms.
In practice, this is further abstracted and users should use the provided term
objects described below or subclass Term
to implement a custom
kernel. Subclasses should overload the
Term.get_real_coefficients()
and
Term.get_complex_coefficients()
methods to provide custom
\(\alpha\) and \(\beta\) arrays.
The following gives a simple example for a custom term of the form:
with the parameters log_a
, log_b
, log_c
, and log_P
.
# If we use the autograd wrappers around numpy then the gradients
# of the term will be computed automatically.
import autograd.numpy as np
class CustomTerm(terms.Term):
parameter_names = ("log_a", "log_b", "log_c", "log_P")
def get_real_coefficients(self, params):
log_a, log_b, log_c, log_P = params
b = np.exp(log_b)
return (
np.exp(log_a) * (1.0 + b) / (2.0 + b), np.exp(log_c),
)
def get_complex_coefficients(self, params):
log_a, log_b, log_c, log_P = params
b = np.exp(log_b)
return (
np.exp(log_a) / (2.0 + b), 0.0,
np.exp(log_c), 2*np.pi*np.exp(-log_P),
)
In this example, all of the returned coefficients are scalars but they can also be returned as arrays. To use this kernel, we would use something like the following:
bounds = dict(log_a=(None, None), log_b=(None, 5.0), log_c=(-1.0, 1.0),
log_P=(-0.5, 0.5))
kernel = CustomTerm(log_a=0.1, log_b=0.5, log_c=-0.01, log_P=0.0,
bounds=bounds)
where bounds
is optional.
Finally, terms are combined by adding or multiplying them:
term1 = terms.RealTerm(0.5, -5.0)
term2 = terms.ComplexTerm(1.0, 0.1, 0.6, -0.3)
term3 = terms.SHOTerm(1.0, 0.1, 0.6)
# Adding kernels works
kernel = term1 + term2
# So does multiplying them
kernel = term1 * term2
# Or, any combination of these operations
kernel = term1 * (term2 + term3)
These combined terms are also Term
objects and they provide the same
interface, concatenating the parameters and coefficients in the correct order.
Since Term
objects implement the Modeling Protocol, for the kernel
described above, we can do things like the following:
print(kernel.get_parameter_dict())
With the output:
OrderedDict([('terms[0]:log_a', 0.5),
('terms[0]:log_c', -5.0),
('terms[1]:log_a', 1.0),
('terms[1]:log_b', 0.1),
('terms[1]:log_c', 0.6),
('terms[1]:log_d', -0.2)])
The following are the Term
objects provided by celerite.
-
class
celerite.terms.
Term
(*args, **kwargs)¶ The abstract base “term” that is the superclass of all other terms
Subclasses should overload the
terms.Term.get_real_coefficients()
andterms.Term.get_complex_coefficients()
methods.-
check_parameters
()¶ Check for negative power in the PSD using Sturm’s theorem
- Returns
True
for valid parameters.
-
property
coefficients
¶ All of the coefficient arrays
This property is the concatenation of the results from
terms.Term.get_real_coefficients()
andterms.Term.get_complex_coefficients()
but it will always return a tuple of length 6, even ifalpha_complex_imag
was omitted fromget_complex_coefficients
.- Returns
(array[j_real], array[j_real], array[j_complex], array[j_complex], array[j_complex], array[j_complex]):
alpha_real
,beta_real
,alpha_complex_real
,alpha_complex_imag
,beta_complex_real
, andbeta_complex_imag
as described above.- Raises
ValueError – For invalid dimensions for the coefficients.
-
get_complex_coefficients
(params)¶ Get the arrays
alpha_complex_*
andbeta_complex_*
This method should be overloaded by subclasses to return the arrays
alpha_complex_real
,alpha_complex_imag
,beta_complex_real
, andbeta_complex_imag
given the current parameter settings. By default, this term is empty.- Returns
(array[j_complex], array[j_complex], array[j_complex], array[j_complex]):
alpha_complex_real
,alpha_complex_imag
,beta_complex_real
, andbeta_complex_imag
as described above.alpha_complex_imag
can be omitted and it will be assumed to be zero.
-
get_psd
(omega)¶ Compute the PSD of the term for an array of angular frequencies
- Parameters
omega (array[..]) – An array of frequencies where the PSD should be evaluated.
- Returns
The value of the PSD for each
omega
. This will have the same shape asomega
.
-
get_real_coefficients
(params)¶ Get the arrays
alpha_real
andbeta_real
This method should be overloaded by subclasses to return the arrays
alpha_real
andbeta_real
given the current parameter settings. By default, this term is empty.- Returns
alpha_real
andbeta_real
as described above.- Return type
(array[j_real], array[j_real])
-
get_value
(tau)¶ Compute the value of the term for an array of lags
- Parameters
tau (array[..]) – An array of lags where the term should be evaluated.
- Returns
The value of the term for each
tau
. This will have the same shape astau
.
-
property
terms
¶ A list of all the terms included in a sum of terms
-
-
class
celerite.terms.
RealTerm
(*args, **kwargs)¶ The simplest celerite term
This term has the form
\[k(\tau) = a_j\,e^{-c_j\,\tau}\]with the parameters
log_a
andlog_c
.Strictly speaking, for a sum of terms, the parameter
a
could be allowed to go negative but since it is somewhat subtle to ensure positive definiteness, we require that the amplitude be positive through this interface. Advanced users can build a custom term that has negative coefficients but care should be taken to ensure positivity.
-
class
celerite.terms.
ComplexTerm
(*args, **kwargs)¶ A general celerite term
This term has the form
\[k(\tau) = \frac{1}{2}\,\left[(a_j + b_j)\,e^{-(c_j+d_j)\,\tau} + (a_j - b_j)\,e^{-(c_j-d_j)\,\tau}\right]\]with the parameters
log_a
,log_b
,log_c
, andlog_d
. The parameterlog_b
can be omitted and it will be assumed to be zero.This term will only correspond to a positive definite kernel (on its own) if \(a_j\,c_j \ge b_j\,d_j\) and the
log_prior
method checks for this constraint.
-
class
celerite.terms.
SHOTerm
(*args, **kwargs)¶ A term representing a stochastically-driven, damped harmonic oscillator
The PSD of this term is
\[S(\omega) = \sqrt{\frac{2}{\pi}} \frac{S_0\,\omega_0^4} {(\omega^2-{\omega_0}^2)^2 + {\omega_0}^2\,\omega^2/Q^2}\]with the parameters
log_S0
,log_Q
, andlog_omega0
.
-
class
celerite.terms.
Matern32Term
(*args, **kwargs)¶ A term that approximates a Matern-3/2 function
This term is defined as
\[k(\tau) = \sigma^2\,\left[ \left(1+1/\epsilon\right)\,e^{-(1-\epsilon)\sqrt{3}\,\tau/\rho} \left(1-1/\epsilon\right)\,e^{-(1+\epsilon)\sqrt{3}\,\tau/\rho} \right]\]with the parameters
log_sigma
andlog_rho
. The parametereps
controls the quality of the approximation since, in the limit \(\epsilon \to 0\) this becomes the Matern-3/2 function\[\lim_{\epsilon \to 0} k(\tau) = \sigma^2\,\left(1+ \frac{\sqrt{3}\,\tau}{\rho}\right)\, \exp\left(-\frac{\sqrt{3}\,\tau}{\rho}\right)\]
-
class
celerite.terms.
JitterTerm
(*args, **kwargs)¶ A diagonal jitter or “white noise” term
This term has the form
\[k(\tau_{n,m}) = \sigma^2\,\delta_{n,m}\]with the parameter
log_sigma
.Note
Jitter is never used in
celerite.GP.predict()
and orcelerite.terms.Term.get_psd()
. If you want to compute the jitter power, it isterm.jitter * N
, whereN
is the number of data points.- Parameters
log_sigma (float) – The log of the amplitude of the white noise.