bcam.indar.ema.LTIKernel¶

class bcam.indar.ema.LTIKernel(*, dt: float = 1.0, alpha: float = 0.0, beta: float = 0.0, mode: str = 'g', atol: float = 1e-06, btol: float = 1e-06, maxiter: int = None, conlim: float = 100000000.0, show: bool = False)¶

Fit a Linear Time-Invariant (LTI) kernel.

In a discrete LTI model, given an input \(x\), the response \(y\) at a given time \(k\) is

\[y_k = dt\,\sum_{i=0}^k h_{k-i}\,x_i + \varepsilon_k,\]

where \(dt\) is the time step, \(h\) is the Impulse Response Function (IRF) to be estimated (also known as the kernel), and \(\varepsilon\) is white noise. This class only admits SISO data, that is, the input and response are scalar time series, and the kernel is a 1d-array.

This class uses the matrix-free algorithm LSMR, as implemented in scipy.sparse.linalg.lsmr(), to estimate the IRF.

Parameters:

dtfloat, default=1.0: Sampling time step.
alpha, betafloat, default=0.0: If alpha is positive, the \(\ell^2\) norm of the kernel is penalized. If beta is positive, the \(\ell^2\) norm of the kernel’s derivative is penalized.
mode{‘g’, ‘a’}, default=’g’: When the mode is general (g), the derivative is penalized for all times, while for acceleration (a) mode, time \(k = 0\) is not penalized.
atol, btolfloat, default=1e-6: atol is the relative tolerance in the entries of the input \(x\), and btol is the relative tolerance in the entries of the response \(y\).
maxiter, conlim, showint, float, bool: See scipy.sparse.linalg.lsmr() for details.

Attributes:

kernel_np.ndarray of shape (n_samples,)

Estimated IRF, or kernel.

info_dict

Information about the optimization process, containing the following keys:

istop: reason for stopping.
itn: number of iterations.
normr: the norm of the residual.
normar: the norm of the projected residual.
norma: the estimate of the Frobenius norm of the matrix.
conda: the estimate of the condition number of the matrix.
normx: the norm of the solution.

See scipy.sparse.linalg.lsmr() for details.

Methods

`fit`(X, y)	Fit LTI model.
`get_metadata_routing`()	Get metadata routing of this object.
`get_params`([deep])	Get parameters for this estimator.
`predict`(X)	Predict response using the LTI model.
`score`(X, y)	Return coefficient of determination on test data.
`set_params`(**params)	Set the parameters of this estimator.

Notes

This class is a sklearn predictor.

Examples

We create synthetic data from a known kernel and fit the LTI model to check if the kernel is correctly estimated.

>>> import numpy as np
>>> from scipy.signal import fftconvolve
>>> from bcam.indar import ema
>>> # Create random input data.
>>> rng = np.random.default_rng(123)
>>> n_reps = 5 # Number of repetitions (trials)
>>> N = 100 # Number of time samples
>>> X = rng.normal(0, 1, size=(n_reps, N))
>>> # Create a kernel and generate response data.
>>> true_kernel = np.exp(-0.05*np.arange(N)) * np.cos(2*np.pi*0.1*np.arange(N))
>>> y = fftconvolve(
...     X, true_kernel[np.newaxis, :], mode='full', axes=1)[:, :N]
>>> # Add noise to the response.
>>> y += rng.normal(0, .1, y.shape)
>>> # Fit kernel from input/output data.
>>> model = ema.LTIKernel(beta=5.)
>>> model.fit(X, y)

We check that the predicted response matches the true response.

>>> import matplotlib.pyplot as plt
>>> plt.plot(true_kernel, label='True kernel')
>>> plt.plot(model.kernel_, label='Estimated kernel')
>>> plt.legend()
>>> plt.show()

(png)

../_images/bcam-indar-ema-LTIKernel-2.png

fit(X: array - like, y: array - like) → LTIKernel¶

Fit LTI model.

Parameters:

Xarray-like of shape (n_repetitions, n_time_samples): Input data.
yarray-like of shape (n_repetitions, n_time_samples): Response data with the same shape as X.

Returns:

selfLTIKernel: Fitted estimator.

get_metadata_routing()¶

Get metadata routing of this object.

Please check User Guide on how the routing mechanism works.

Returns:

routingMetadataRequest: A MetadataRequest encapsulating routing information.

get_params(deep=True)¶

Get parameters for this estimator.

Parameters:

deepbool, default=True: If True, will return the parameters for this estimator and contained subobjects that are estimators.

Returns:

paramsdict: Parameter names mapped to their values.

predict(X: array - like) → ndarray¶

Predict response using the LTI model.

Parameters:

Xarray-like of shape (n_repetitions, n_time_samples): Input data.

Returns:

y_prednp.ndarray of shape (n_repetitions, n_time_samples): Predicted response data.

score(X: array - like, y: array - like) → float¶

Return coefficient of determination on test data.

Parameters:

Xarray-like of shape (n_repetitions, n_time_samples): Input data.
yarray-like of shape (n_repetitions, n_time_samples): True response data.

Returns:

scorefloat: \(R^2\) of self.predict(X) with respect to y.

set_params(**params)¶

Set the parameters of this estimator.

The method works on simple estimators as well as on nested objects (such as Pipeline). The latter have parameters of the form <component>__<parameter> so that it’s possible to update each component of a nested object.

Parameters:

**paramsdict: Estimator parameters.

Returns:

selfestimator instance: Estimator instance.