Source code for gammapy.modeling.covariance

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Covariance class"""
import numpy as np
import scipy
from .parameter import Parameters

__all__ = ["Covariance"]


[docs]class Covariance: """Parameter covariance class Parameters ---------- parameters : `~gammapy.modeling.Parameters` Parameter list data : `~numpy.ndarray` Covariance data array """ def __init__(self, parameters, data=None): self.parameters = parameters if data is None: data = np.diag([p.error ** 2 for p in self.parameters]) self._data = np.asanyarray(data, dtype=float) @property def shape(self): """Covariance shape""" npars = len(self.parameters) return npars, npars @property def data(self): """Covariance data (`~numpy.ndarray`)""" return self._data @data.setter def data(self, value): value = np.asanyarray(value) npars = len(self.parameters) shape = (npars, npars) if value.shape != shape: raise ValueError( f"Invalid covariance shape: {value.shape}, expected {shape}" ) self._data = value @staticmethod def _expand_factor_matrix(matrix, parameters): """Expand covariance matrix with zeros for frozen parameters""" npars = len(parameters) matrix_expanded = np.zeros((npars, npars)) mask_frozen = [par.frozen for par in parameters] pars_index = [np.where(np.array(parameters) == p)[0][0] for p in parameters] mask_duplicate = [pars_idx != idx for idx, pars_idx in enumerate(pars_index)] mask = np.array(mask_frozen) | np.array(mask_duplicate) free_parameters = ~(mask | mask[:, np.newaxis]) matrix_expanded[free_parameters] = matrix.ravel() return matrix_expanded
[docs] @classmethod def from_factor_matrix(cls, parameters, matrix): """Set covariance from factor covariance matrix. Used in the optimizer interface. """ npars = len(parameters) if not matrix.shape == (npars, npars): matrix = cls._expand_factor_matrix(matrix, parameters) scales = [par.scale for par in parameters] scale_matrix = np.outer(scales, scales) data = scale_matrix * matrix return cls(parameters, data=data)
[docs] @classmethod def from_stack(cls, covar_list): """Stack sub-covariance matrices from list Parameters ---------- covar_list : list of `Covariance` List of sub-covariances Returns ------- covar : `Covariance` Stacked covariance """ parameters = Parameters.from_stack([_.parameters for _ in covar_list]) covar = cls(parameters) for subcovar in covar_list: covar.set_subcovariance(subcovar) return covar
[docs] def get_subcovariance(self, parameters): """Get sub-covariance matrix Parameters ---------- parameters : `Parameters` Sub list of parameters. Returns ------- covariance : `~numpy.ndarray` Sub-covariance. """ idx = [self.parameters.index(par) for par in parameters] data = self._data[np.ix_(idx, idx)] return self.__class__(parameters=parameters, data=data)
[docs] def set_subcovariance(self, covar): """Set sub-covariance matrix Parameters ---------- parameters : `Parameters` Sub list of parameters. """ idx = [self.parameters.index(par) for par in covar.parameters] if not np.allclose(self.data[np.ix_(idx, idx)], covar.data): self.data[idx, :] = 0 self.data[:, idx] = 0 self._data[np.ix_(idx, idx)] = covar.data
[docs] def plot_correlation(self, ax=None, **kwargs): """Plot correlation matrix. Parameters ---------- ax : `~matplotlib.axes.Axes`, optional Axis to plot on. **kwargs : dict Keyword arguments passed to `~gammapy.visualisation.plot_heatmap` Returns ------- ax : `~matplotlib.axes.Axes`, optional Axis """ import matplotlib.pyplot as plt from gammapy.visualization import plot_heatmap, annotate_heatmap npars = len(self.parameters) figsize = (npars * 0.8, npars * 0.65) plt.figure(figsize=figsize) ax = plt.gca() if ax is None else ax kwargs.setdefault("cmap", "coolwarm") names = self.parameters.names im, cbar = plot_heatmap( data=self.correlation, col_labels=names, row_labels=names, ax=ax, vmin=-1, vmax=1, cbarlabel="Correlation", **kwargs, ) annotate_heatmap(im=im) return ax
@property def correlation(self): r"""Correlation matrix (`numpy.ndarray`). Correlation :math:`C` is related to covariance :math:`\Sigma` via: .. math:: C_{ij} = \frac{ \Sigma_{ij} }{ \sqrt{\Sigma_{ii} \Sigma_{jj}} } """ err = np.sqrt(np.diag(self.data)) with np.errstate(invalid="ignore", divide="ignore"): correlation = self.data / np.outer(err, err) return np.nan_to_num(correlation) @property def scipy_mvn(self): # TODO: use this, as in https://github.com/cdeil/multinorm/blob/master/multinorm.py return scipy.stats.multivariate_normal( self.parameters.values, self.data, allow_singular=True ) def __str__(self): return str(self.data) def __array__(self): return self.data