Module gmr.mvn
Functions
def condition(mean, covariance, i_out, i_in, X)
-
Expand source code
def condition(mean, covariance, i_out, i_in, X): """Compute conditional mean and covariance. Parameters ---------- mean : array, shape (n_features,) Mean of MVN covariance : array, shape (n_features, n_features) Covariance of MVN i_out : array, shape (n_features_out,) Output feature indices i_in : array, shape (n_features_in,) Input feature indices X : array, shape (n_samples, n_features_out) Inputs Returns ------- mean : array, shape (n_features_out,) Mean of the conditional distribution covariance : array, shape (n_features_out, n_features_out) Covariance of the conditional distribution """ cov_12 = covariance[np.ix_(i_out, i_in)] cov_11 = covariance[np.ix_(i_out, i_out)] regression_coeffs = regression_coefficients( covariance, i_out, i_in, cov_12=cov_12) mean = mean[i_out] + regression_coeffs.dot((X - mean[i_in]).T).T covariance = cov_11 - regression_coeffs.dot(cov_12.T) return mean, covariance
Compute conditional mean and covariance.
Parameters
mean
:array, shape (n_features,)
- Mean of MVN
covariance
:array, shape (n_features, n_features)
- Covariance of MVN
i_out
:array, shape (n_features_out,)
- Output feature indices
i_in
:array, shape (n_features_in,)
- Input feature indices
X
:array, shape (n_samples, n_features_out)
- Inputs
Returns
mean
:array, shape (n_features_out,)
- Mean of the conditional distribution
covariance
:array, shape (n_features_out, n_features_out)
- Covariance of the conditional distribution
def invert_indices(n_features, indices)
-
Expand source code
def invert_indices(n_features, indices): inv = np.ones(n_features, dtype=bool) inv[indices] = False inv, = np.where(inv) return inv
def plot_error_ellipse(ax,
mvn,
color=None,
alpha=0.25,
factors=array([0.25, 0.5 , 0.75, 1. , 1.25, 1.5 , 1.75, 2. ]))-
Expand source code
def plot_error_ellipse(ax, mvn, color=None, alpha=0.25, factors=np.linspace(0.25, 2.0, 8)): """Plot error ellipse of MVN. Parameters ---------- ax : axis Matplotlib axis. mvn : MVN Multivariate normal distribution. color : str, optional (default: None) Color in which the ellipse should be plotted alpha : float, optional (default: 0.25) Alpha value for ellipse factors : array, optional (default: np.linspace(0.25, 2.0, 8)) Multiples of the standard deviations that should be plotted. """ from matplotlib.patches import Ellipse for factor in factors: angle, width, height = mvn.to_ellipse(factor) ell = Ellipse(xy=mvn.mean, width=2.0 * width, height=2.0 * height, angle=np.degrees(angle)) ell.set_alpha(alpha) if color is not None: ell.set_color(color) ax.add_artist(ell)
Plot error ellipse of MVN.
Parameters
ax
:axis
- Matplotlib axis.
mvn
:MVN
- Multivariate normal distribution.
color
:str
, optional(default: None)
- Color in which the ellipse should be plotted
alpha
:float
, optional(default: 0.25)
- Alpha value for ellipse
factors
:array
, optional(default: np.linspace(0.25, 2.0, 8))
- Multiples of the standard deviations that should be plotted.
def regression_coefficients(covariance, i_out, i_in, cov_12=None)
-
Expand source code
def regression_coefficients(covariance, i_out, i_in, cov_12=None): """Compute regression coefficients to predict conditional distribution. Parameters ---------- covariance : array, shape (n_features, n_features) Covariance of MVN i_out : array, shape (n_features_out,) Output feature indices i_in : array, shape (n_features_in,) Input feature indices cov_12 : array, shape (n_features_out, n_features_in), optional (default: None) Precomputed block of the covariance matrix between input features and output features Returns ------- regression_coeffs : array, shape (n_features_out, n_features_in) Regression coefficients. These can be used to compute the mean of the conditional distribution as mean[i1] + regression_coeffs.dot((X - mean[i2]).T).T """ if cov_12 is None: cov_12 = covariance[np.ix_(i_out, i_in)] cov_22 = covariance[np.ix_(i_in, i_in)] prec_22 = pinvh(cov_22) return cov_12.dot(prec_22)
Compute regression coefficients to predict conditional distribution.
Parameters
covariance
:array, shape (n_features, n_features)
- Covariance of MVN
i_out
:array, shape (n_features_out,)
- Output feature indices
i_in
:array, shape (n_features_in,)
- Input feature indices
cov_12
:array, shape (n_features_out, n_features_in)
, optional(default: None)
- Precomputed block of the covariance matrix between input features and output features
Returns
regression_coeffs
:array, shape (n_features_out, n_features_in)
- Regression coefficients. These can be used to compute the mean of the conditional distribution as mean[i1] + regression_coeffs.dot((X - mean[i2]).T).T
Classes
class MVN (mean=None, covariance=None, verbose=0, random_state=None)
-
Expand source code
class MVN(object): """Multivariate normal distribution. Some utility functions for MVNs. See http://en.wikipedia.org/wiki/Multivariate_normal_distribution for more details. Parameters ---------- mean : array-like, shape (n_features), optional Mean of the MVN. covariance : array-like, shape (n_features, n_features), optional Covariance of the MVN. verbose : int, optional (default: 0) Verbosity level. random_state : int or RandomState, optional (default: global random state) If an integer is given, it fixes the seed. Defaults to the global numpy random number generator. """ def __init__(self, mean=None, covariance=None, verbose=0, random_state=None): self.mean = mean self.covariance = covariance self.verbose = verbose self.random_state = check_random_state(random_state) self.norm = None if self.mean is not None: self.mean = np.asarray(self.mean) if self.covariance is not None: self.covariance = np.asarray(self.covariance) def _check_initialized(self): if self.mean is None: raise ValueError("Mean has not been initialized") if self.covariance is None: raise ValueError("Covariance has not been initialized") def from_samples(self, X, bessels_correction=True): """MLE of the mean and covariance. Parameters ---------- X : array-like, shape (n_samples, n_features) Samples from the true function. bessels_correction : bool Apply Bessel's correction to the covariance estimate. Returns ------- self : MVN This object. """ self.mean = np.mean(X, axis=0) bias = 0 if bessels_correction else 1 self.covariance = np.cov(X, rowvar=0, bias=bias) self.norm = None return self def sample(self, n_samples): """Sample from multivariate normal distribution. Parameters ---------- n_samples : int Number of samples. Returns ------- X : array, shape (n_samples, n_features) Samples from the MVN. """ self._check_initialized() return self.random_state.multivariate_normal( self.mean, self.covariance, size=(n_samples,)) def sample_confidence_region(self, n_samples, alpha): """Sample from alpha confidence region. Parameters ---------- n_samples : int Number of samples. alpha : float Value between 0 and 1 that defines the probability of the confidence region, e.g., 0.6827 for the 1-sigma confidence region or 0.9545 for the 2-sigma confidence region. Returns ------- X : array, shape (n_samples, n_features) Samples from the confidence region. """ return np.array([self._one_sample_confidence_region(alpha) for _ in range(n_samples)]) def _one_sample_confidence_region(self, alpha): x = self.sample(1)[0] while not self.is_in_confidence_region(x, alpha): x = self.sample(1)[0] return x def is_in_confidence_region(self, x, alpha): """Check if sample is in alpha confidence region. Parameters ---------- x : array, shape (n_features,) Sample alpha : float Value between 0 and 1 that defines the probability of the confidence region, e.g., 0.6827 for the 1-sigma confidence region or 0.9545 for the 2-sigma confidence region. Returns ------- is_in_confidence_region : bool Is the sample in the alpha confidence region? """ self._check_initialized() # we have one degree of freedom less than number of dimensions n_dof = len(x) - 1 if n_dof >= 1: return self.squared_mahalanobis_distance(x) <= chi2(n_dof).ppf(alpha) else: # 1D lo, hi = norm.interval( alpha, loc=self.mean[0], scale=self.covariance[0, 0]) return lo <= x[0] <= hi def to_norm_factor_and_exponents(self, X): """Compute normalization factor and exponents of Gaussian. These values can be used to compute the probability density function of this Gaussian: p(x) = norm_factor * np.exp(exponents). Parameters ---------- X : array-like, shape (n_samples, n_features) Data. Returns ------- norm_factor : float Normalization factor: constant term outside of exponential function in probability density function of this Gaussian. exponents : array, shape (n_samples,) Exponents to compute probability density function. """ self._check_initialized() X = np.atleast_2d(X) n_features = X.shape[1] try: L = sp.linalg.cholesky(self.covariance, lower=True) except np.linalg.LinAlgError: # Degenerated covariance, try to add regularization L = sp.linalg.cholesky( self.covariance + 1e-3 * np.eye(n_features), lower=True) X_minus_mean = X - self.mean if self.norm is None: # Suppress a determinant of 0 to avoid numerical problems L_det = max(sp.linalg.det(L), np.finfo(L.dtype).eps) self.norm = (2.0 * np.pi) ** (-0.5 * n_features) / L_det # Solve L x = (X - mean)^T for x with triangular L # (LL^T = Sigma), that is, x = L^T^-1 (X - mean)^T. # We can avoid covariance inversion when computing # (X - mean) Sigma^-1 (X - mean)^T with this trick, # since Sigma^-1 = L^T^-1 L^-1. X_normalized = sp.linalg.solve_triangular( L, X_minus_mean.T, lower=True).T exponent = -0.5 * np.sum(X_normalized ** 2, axis=1) return self.norm, exponent def to_probability_density(self, X): """Compute probability density. Parameters ---------- X : array-like, shape (n_samples, n_features) Data. Returns ------- p : array, shape (n_samples,) Probability densities of data. """ norm_factor, exponents = self.to_norm_factor_and_exponents(X) return norm_factor * np.exp(exponents) def marginalize(self, indices): """Marginalize over everything except the given indices. Parameters ---------- indices : array, shape (n_new_features,) Indices of dimensions that we want to keep. Returns ------- marginal : MVN Marginal MVN distribution. """ self._check_initialized() return MVN(mean=self.mean[indices], covariance=self.covariance[np.ix_(indices, indices)]) def condition(self, indices, x): """Conditional distribution over given indices. Parameters ---------- indices : array, shape (n_new_features,) Indices of dimensions that we want to condition. x : array, shape (n_new_features,) Values of the features that we know. Returns ------- conditional : MVN Conditional MVN distribution p(Y | X=x). """ self._check_initialized() mean, covariance = condition( self.mean, self.covariance, invert_indices(self.mean.shape[0], indices), indices, x) return MVN(mean=mean, covariance=covariance, random_state=self.random_state) def predict(self, indices, X): """Predict means and covariance of posteriors. Same as condition() but for multiple samples. Parameters ---------- indices : array-like, shape (n_features_in,) Indices of dimensions that we want to condition. X : array-like, shape (n_samples, n_features_in) Values of the features that we know. Returns ------- Y : array, shape (n_samples, n_features_out) Predicted means of missing values. covariance : array, shape (n_features_out, n_features_out) Covariance of the predicted features. """ self._check_initialized() indices = np.asarray(indices, dtype=int) X = np.asarray(X) return condition( self.mean, self.covariance, invert_indices(self.mean.shape[0], indices), indices, X) def squared_mahalanobis_distance(self, x): """Squared Mahalanobis distance between point and this MVN. Parameters ---------- x : array, shape (n_features,) Returns ------- d : float Squared Mahalanobis distance """ self._check_initialized() return mahalanobis(x, self.mean, np.linalg.inv(self.covariance)) ** 2 def to_ellipse(self, factor=1.0): """Compute error ellipse. An error ellipse shows equiprobable points. Parameters ---------- factor : float One means standard deviation. Returns ------- angle : float Rotation angle of the ellipse. width : float Width of the ellipse (semi axis, not diameter). height : float Height of the ellipse (semi axis, not diameter). """ self._check_initialized() vals, vecs = sp.linalg.eigh(self.covariance) order = vals.argsort()[::-1] vals, vecs = vals[order], vecs[:, order] angle = np.arctan2(*vecs[:, 0][::-1]) width, height = factor * np.sqrt(vals) return angle, width, height def _sqrt_cov(self, C): """Compute square root of a symmetric matrix. Parameters ---------- C : array, shape (n_features, n_features) Symmetric matrix. Returns ------- sqrt(C) : array, shape (n_features, n_features) Square root of covariance. The square root of a square matrix is defined as :math:`\Sigma^{\frac{1}{2}} = B \sqrt(D) B^T`, where :math:`\Sigma = B D B^T` is the Eigen decomposition of the covariance. """ D, B = np.linalg.eigh(C) # HACK: avoid numerical problems D = np.maximum(D, np.finfo(float).eps) return B.dot(np.diag(np.sqrt(D))).dot(B.T) def sigma_points(self, alpha=1e-3, kappa=0.0): """Compute sigma points for unscented transform. The unscented transform allows us to estimate the resulting MVN from applying a nonlinear transformation :math:`f` to this MVN. In order to do this, you have to transform the sigma points obtained from this function with :math:`f` and then create the new MVN with :func:`MVN.estimate_from_sigma_points`. The unscented transform is most commonly used in the Unscented Kalman Filter (UKF). Parameters ---------- alpha : float, optional (default: 1e-3) Determines the spread of the sigma points around the mean and is usually set to a small positive value. kappa : float, optional (default: 0) A secondary scaling parameter which is usually set to 0. Returns ------- sigma_points : array, shape (2 * n_features + 1, n_features) Query points that have to be transformed to estimate the resulting MVN. """ self._check_initialized() n_features = len(self.mean) lmbda = alpha ** 2 * (n_features + kappa) - n_features offset = self._sqrt_cov((n_features + lmbda) * self.covariance) points = np.empty(((2 * n_features + 1), n_features)) points[0, :] = self.mean for i in range(n_features): points[1 + i, :] = self.mean + offset[i] points[1 + n_features + i:, :] = self.mean - offset[i] return points def estimate_from_sigma_points(self, transformed_sigma_points, alpha=1e-3, beta=2.0, kappa=0.0, random_state=None): """Estimate new MVN from sigma points through the unscented transform. See :func:`MVN.sigma_points` for more details. Parameters ---------- transformed_sigma_points : array, shape (2 * n_features + 1, n_features) Query points that were transformed to estimate the resulting MVN. alpha : float, optional (default: 1e-3) Determines the spread of the sigma points around the mean and is usually set to a small positive value. Note that this value has to match the value that was used to create the sigma points. beta : float, optional (default: 2) Encodes information about the distribution. For Gaussian distributions, beta=2 is the optimal choice. kappa : float, optional (default: 0) A secondary scaling parameter which is usually set to 0. Note that this value has to match the value that was used to create the sigma points. random_state : int or RandomState, optional (default: random state of self) If an integer is given, it fixes the seed. Defaults to the global numpy random number generator. Returns ------- mvn : MVN Transformed MVN: f(self). """ self._check_initialized() n_features = len(self.mean) lmbda = alpha ** 2 * (n_features + kappa) - n_features mean_weight_0 = lmbda / (n_features + lmbda) cov_weight_0 = lmbda / (n_features + lmbda) + (1 - alpha ** 2 + beta) weights_i = 1.0 / (2.0 * (n_features + lmbda)) mean_weights = np.empty(len(transformed_sigma_points)) mean_weights[0] = mean_weight_0 mean_weights[1:] = weights_i cov_weights = np.empty(len(transformed_sigma_points)) cov_weights[0] = cov_weight_0 cov_weights[1:] = weights_i mean = np.sum(mean_weights[:, np.newaxis] * transformed_sigma_points, axis=0) sigma_points_minus_mean = transformed_sigma_points - mean covariance = sigma_points_minus_mean.T.dot( np.diag(cov_weights)).dot(sigma_points_minus_mean) if random_state is None: random_state = self.random_state return MVN(mean=mean, covariance=covariance, random_state=random_state)
Multivariate normal distribution.
Some utility functions for MVNs. See http://en.wikipedia.org/wiki/Multivariate_normal_distribution for more details.
Parameters
mean
:array-like, shape (n_features)
, optional- Mean of the MVN.
covariance
:array-like, shape (n_features, n_features)
, optional- Covariance of the MVN.
verbose
:int
, optional(default: 0)
- Verbosity level.
random_state
:int
orRandomState
, optional(default: global random state)
- If an integer is given, it fixes the seed. Defaults to the global numpy random number generator.
Methods
def condition(self, indices, x)
-
Expand source code
def condition(self, indices, x): """Conditional distribution over given indices. Parameters ---------- indices : array, shape (n_new_features,) Indices of dimensions that we want to condition. x : array, shape (n_new_features,) Values of the features that we know. Returns ------- conditional : MVN Conditional MVN distribution p(Y | X=x). """ self._check_initialized() mean, covariance = condition( self.mean, self.covariance, invert_indices(self.mean.shape[0], indices), indices, x) return MVN(mean=mean, covariance=covariance, random_state=self.random_state)
Conditional distribution over given indices.
Parameters
indices
:array, shape (n_new_features,)
- Indices of dimensions that we want to condition.
x
:array, shape (n_new_features,)
- Values of the features that we know.
Returns
conditional
:MVN
- Conditional MVN distribution p(Y | X=x).
def estimate_from_sigma_points(self, transformed_sigma_points, alpha=0.001, beta=2.0, kappa=0.0, random_state=None)
-
Expand source code
def estimate_from_sigma_points(self, transformed_sigma_points, alpha=1e-3, beta=2.0, kappa=0.0, random_state=None): """Estimate new MVN from sigma points through the unscented transform. See :func:`MVN.sigma_points` for more details. Parameters ---------- transformed_sigma_points : array, shape (2 * n_features + 1, n_features) Query points that were transformed to estimate the resulting MVN. alpha : float, optional (default: 1e-3) Determines the spread of the sigma points around the mean and is usually set to a small positive value. Note that this value has to match the value that was used to create the sigma points. beta : float, optional (default: 2) Encodes information about the distribution. For Gaussian distributions, beta=2 is the optimal choice. kappa : float, optional (default: 0) A secondary scaling parameter which is usually set to 0. Note that this value has to match the value that was used to create the sigma points. random_state : int or RandomState, optional (default: random state of self) If an integer is given, it fixes the seed. Defaults to the global numpy random number generator. Returns ------- mvn : MVN Transformed MVN: f(self). """ self._check_initialized() n_features = len(self.mean) lmbda = alpha ** 2 * (n_features + kappa) - n_features mean_weight_0 = lmbda / (n_features + lmbda) cov_weight_0 = lmbda / (n_features + lmbda) + (1 - alpha ** 2 + beta) weights_i = 1.0 / (2.0 * (n_features + lmbda)) mean_weights = np.empty(len(transformed_sigma_points)) mean_weights[0] = mean_weight_0 mean_weights[1:] = weights_i cov_weights = np.empty(len(transformed_sigma_points)) cov_weights[0] = cov_weight_0 cov_weights[1:] = weights_i mean = np.sum(mean_weights[:, np.newaxis] * transformed_sigma_points, axis=0) sigma_points_minus_mean = transformed_sigma_points - mean covariance = sigma_points_minus_mean.T.dot( np.diag(cov_weights)).dot(sigma_points_minus_mean) if random_state is None: random_state = self.random_state return MVN(mean=mean, covariance=covariance, random_state=random_state)
Estimate new MVN from sigma points through the unscented transform.
See :func:
MVN.sigma_points()
for more details.Parameters
transformed_sigma_points
:array, shape (2 * n_features + 1, n_features)
- Query points that were transformed to estimate the resulting MVN.
alpha
:float
, optional(default: 1e-3)
- Determines the spread of the sigma points around the mean and is usually set to a small positive value. Note that this value has to match the value that was used to create the sigma points.
beta
:float
, optional(default: 2)
- Encodes information about the distribution. For Gaussian distributions, beta=2 is the optimal choice.
kappa
:float
, optional(default: 0)
- A secondary scaling parameter which is usually set to 0. Note that this value has to match the value that was used to create the sigma points.
random_state
:int
orRandomState
, optional(default: random state
ofself)
- If an integer is given, it fixes the seed. Defaults to the global numpy random number generator.
Returns
mvn
:MVN
- Transformed MVN: f(self).
def from_samples(self, X, bessels_correction=True)
-
Expand source code
def from_samples(self, X, bessels_correction=True): """MLE of the mean and covariance. Parameters ---------- X : array-like, shape (n_samples, n_features) Samples from the true function. bessels_correction : bool Apply Bessel's correction to the covariance estimate. Returns ------- self : MVN This object. """ self.mean = np.mean(X, axis=0) bias = 0 if bessels_correction else 1 self.covariance = np.cov(X, rowvar=0, bias=bias) self.norm = None return self
MLE of the mean and covariance.
Parameters
X
:array-like, shape (n_samples, n_features)
- Samples from the true function.
bessels_correction
:bool
- Apply Bessel's correction to the covariance estimate.
Returns
self
:MVN
- This object.
def is_in_confidence_region(self, x, alpha)
-
Expand source code
def is_in_confidence_region(self, x, alpha): """Check if sample is in alpha confidence region. Parameters ---------- x : array, shape (n_features,) Sample alpha : float Value between 0 and 1 that defines the probability of the confidence region, e.g., 0.6827 for the 1-sigma confidence region or 0.9545 for the 2-sigma confidence region. Returns ------- is_in_confidence_region : bool Is the sample in the alpha confidence region? """ self._check_initialized() # we have one degree of freedom less than number of dimensions n_dof = len(x) - 1 if n_dof >= 1: return self.squared_mahalanobis_distance(x) <= chi2(n_dof).ppf(alpha) else: # 1D lo, hi = norm.interval( alpha, loc=self.mean[0], scale=self.covariance[0, 0]) return lo <= x[0] <= hi
Check if sample is in alpha confidence region.
Parameters
x
:array, shape (n_features,)
- Sample
alpha
:float
- Value between 0 and 1 that defines the probability of the confidence region, e.g., 0.6827 for the 1-sigma confidence region or 0.9545 for the 2-sigma confidence region.
Returns
is_in_confidence_region
:bool
- Is the sample in the alpha confidence region?
def marginalize(self, indices)
-
Expand source code
def marginalize(self, indices): """Marginalize over everything except the given indices. Parameters ---------- indices : array, shape (n_new_features,) Indices of dimensions that we want to keep. Returns ------- marginal : MVN Marginal MVN distribution. """ self._check_initialized() return MVN(mean=self.mean[indices], covariance=self.covariance[np.ix_(indices, indices)])
Marginalize over everything except the given indices.
Parameters
indices
:array, shape (n_new_features,)
- Indices of dimensions that we want to keep.
Returns
marginal
:MVN
- Marginal MVN distribution.
def predict(self, indices, X)
-
Expand source code
def predict(self, indices, X): """Predict means and covariance of posteriors. Same as condition() but for multiple samples. Parameters ---------- indices : array-like, shape (n_features_in,) Indices of dimensions that we want to condition. X : array-like, shape (n_samples, n_features_in) Values of the features that we know. Returns ------- Y : array, shape (n_samples, n_features_out) Predicted means of missing values. covariance : array, shape (n_features_out, n_features_out) Covariance of the predicted features. """ self._check_initialized() indices = np.asarray(indices, dtype=int) X = np.asarray(X) return condition( self.mean, self.covariance, invert_indices(self.mean.shape[0], indices), indices, X)
Predict means and covariance of posteriors.
Same as condition() but for multiple samples.
Parameters
indices
:array-like, shape (n_features_in,)
- Indices of dimensions that we want to condition.
X
:array-like, shape (n_samples, n_features_in)
- Values of the features that we know.
Returns
Y
:array, shape (n_samples, n_features_out)
- Predicted means of missing values.
covariance
:array, shape (n_features_out, n_features_out)
- Covariance of the predicted features.
def sample(self, n_samples)
-
Expand source code
def sample(self, n_samples): """Sample from multivariate normal distribution. Parameters ---------- n_samples : int Number of samples. Returns ------- X : array, shape (n_samples, n_features) Samples from the MVN. """ self._check_initialized() return self.random_state.multivariate_normal( self.mean, self.covariance, size=(n_samples,))
Sample from multivariate normal distribution.
Parameters
n_samples
:int
- Number of samples.
Returns
X
:array, shape (n_samples, n_features)
- Samples from the MVN.
def sample_confidence_region(self, n_samples, alpha)
-
Expand source code
def sample_confidence_region(self, n_samples, alpha): """Sample from alpha confidence region. Parameters ---------- n_samples : int Number of samples. alpha : float Value between 0 and 1 that defines the probability of the confidence region, e.g., 0.6827 for the 1-sigma confidence region or 0.9545 for the 2-sigma confidence region. Returns ------- X : array, shape (n_samples, n_features) Samples from the confidence region. """ return np.array([self._one_sample_confidence_region(alpha) for _ in range(n_samples)])
Sample from alpha confidence region.
Parameters
n_samples
:int
- Number of samples.
alpha
:float
- Value between 0 and 1 that defines the probability of the confidence region, e.g., 0.6827 for the 1-sigma confidence region or 0.9545 for the 2-sigma confidence region.
Returns
X
:array, shape (n_samples, n_features)
- Samples from the confidence region.
def sigma_points(self, alpha=0.001, kappa=0.0)
-
Expand source code
def sigma_points(self, alpha=1e-3, kappa=0.0): """Compute sigma points for unscented transform. The unscented transform allows us to estimate the resulting MVN from applying a nonlinear transformation :math:`f` to this MVN. In order to do this, you have to transform the sigma points obtained from this function with :math:`f` and then create the new MVN with :func:`MVN.estimate_from_sigma_points`. The unscented transform is most commonly used in the Unscented Kalman Filter (UKF). Parameters ---------- alpha : float, optional (default: 1e-3) Determines the spread of the sigma points around the mean and is usually set to a small positive value. kappa : float, optional (default: 0) A secondary scaling parameter which is usually set to 0. Returns ------- sigma_points : array, shape (2 * n_features + 1, n_features) Query points that have to be transformed to estimate the resulting MVN. """ self._check_initialized() n_features = len(self.mean) lmbda = alpha ** 2 * (n_features + kappa) - n_features offset = self._sqrt_cov((n_features + lmbda) * self.covariance) points = np.empty(((2 * n_features + 1), n_features)) points[0, :] = self.mean for i in range(n_features): points[1 + i, :] = self.mean + offset[i] points[1 + n_features + i:, :] = self.mean - offset[i] return points
Compute sigma points for unscented transform.
The unscented transform allows us to estimate the resulting MVN from applying a nonlinear transformation :math:
f
to this MVN. In order to do this, you have to transform the sigma points obtained from this function with :math:f
and then create the new MVN with :func:MVN.estimate_from_sigma_points()
. The unscented transform is most commonly used in the Unscented Kalman Filter (UKF).Parameters
alpha
:float
, optional(default: 1e-3)
- Determines the spread of the sigma points around the mean and is usually set to a small positive value.
kappa
:float
, optional(default: 0)
- A secondary scaling parameter which is usually set to 0.
Returns
sigma_points
:array, shape (2 * n_features + 1, n_features)
- Query points that have to be transformed to estimate the resulting MVN.
def squared_mahalanobis_distance(self, x)
-
Expand source code
def squared_mahalanobis_distance(self, x): """Squared Mahalanobis distance between point and this MVN. Parameters ---------- x : array, shape (n_features,) Returns ------- d : float Squared Mahalanobis distance """ self._check_initialized() return mahalanobis(x, self.mean, np.linalg.inv(self.covariance)) ** 2
Squared Mahalanobis distance between point and this MVN.
Parameters
x
:array, shape (n_features,)
Returns
d
:float
- Squared Mahalanobis distance
def to_ellipse(self, factor=1.0)
-
Expand source code
def to_ellipse(self, factor=1.0): """Compute error ellipse. An error ellipse shows equiprobable points. Parameters ---------- factor : float One means standard deviation. Returns ------- angle : float Rotation angle of the ellipse. width : float Width of the ellipse (semi axis, not diameter). height : float Height of the ellipse (semi axis, not diameter). """ self._check_initialized() vals, vecs = sp.linalg.eigh(self.covariance) order = vals.argsort()[::-1] vals, vecs = vals[order], vecs[:, order] angle = np.arctan2(*vecs[:, 0][::-1]) width, height = factor * np.sqrt(vals) return angle, width, height
Compute error ellipse.
An error ellipse shows equiprobable points.
Parameters
factor
:float
- One means standard deviation.
Returns
angle
:float
- Rotation angle of the ellipse.
width
:float
- Width of the ellipse (semi axis, not diameter).
height
:float
- Height of the ellipse (semi axis, not diameter).
def to_norm_factor_and_exponents(self, X)
-
Expand source code
def to_norm_factor_and_exponents(self, X): """Compute normalization factor and exponents of Gaussian. These values can be used to compute the probability density function of this Gaussian: p(x) = norm_factor * np.exp(exponents). Parameters ---------- X : array-like, shape (n_samples, n_features) Data. Returns ------- norm_factor : float Normalization factor: constant term outside of exponential function in probability density function of this Gaussian. exponents : array, shape (n_samples,) Exponents to compute probability density function. """ self._check_initialized() X = np.atleast_2d(X) n_features = X.shape[1] try: L = sp.linalg.cholesky(self.covariance, lower=True) except np.linalg.LinAlgError: # Degenerated covariance, try to add regularization L = sp.linalg.cholesky( self.covariance + 1e-3 * np.eye(n_features), lower=True) X_minus_mean = X - self.mean if self.norm is None: # Suppress a determinant of 0 to avoid numerical problems L_det = max(sp.linalg.det(L), np.finfo(L.dtype).eps) self.norm = (2.0 * np.pi) ** (-0.5 * n_features) / L_det # Solve L x = (X - mean)^T for x with triangular L # (LL^T = Sigma), that is, x = L^T^-1 (X - mean)^T. # We can avoid covariance inversion when computing # (X - mean) Sigma^-1 (X - mean)^T with this trick, # since Sigma^-1 = L^T^-1 L^-1. X_normalized = sp.linalg.solve_triangular( L, X_minus_mean.T, lower=True).T exponent = -0.5 * np.sum(X_normalized ** 2, axis=1) return self.norm, exponent
Compute normalization factor and exponents of Gaussian.
These values can be used to compute the probability density function of this Gaussian: p(x) = norm_factor * np.exp(exponents).
Parameters
X
:array-like, shape (n_samples, n_features)
- Data.
Returns
norm_factor
:float
- Normalization factor: constant term outside of exponential function in probability density function of this Gaussian.
exponents
:array, shape (n_samples,)
- Exponents to compute probability density function.
def to_probability_density(self, X)
-
Expand source code
def to_probability_density(self, X): """Compute probability density. Parameters ---------- X : array-like, shape (n_samples, n_features) Data. Returns ------- p : array, shape (n_samples,) Probability densities of data. """ norm_factor, exponents = self.to_norm_factor_and_exponents(X) return norm_factor * np.exp(exponents)
Compute probability density.
Parameters
X
:array-like, shape (n_samples, n_features)
- Data.
Returns
p
:array, shape (n_samples,)
- Probability densities of data.