pymc.LKJCholeskyCov#

class pymc.LKJCholeskyCov(name, eta, n, sd_dist, *, compute_corr=True, store_in_trace=True, **kwargs)[source]#

Wrapper class for covariance matrix with LKJ distributed correlations.

This defines a distribution over Cholesky decomposed covariance matrices, such that the underlying correlation matrices follow an LKJ distribution [1] and the standard deviations follow an arbitrary distribution specified by the user.

Parameters:
namestr

The name given to the variable in the model.

etatensor_like of float

The shape parameter (eta > 0) of the LKJ distribution. eta = 1 implies a uniform distribution of the correlation matrices; larger values put more weight on matrices with few correlations.

ntensor_like of int

Dimension of the covariance matrix (n > 1).

sd_distDistribution

A positive scalar or vector distribution for the standard deviations, created with the .dist() API. Should have shape[-1]=n. Scalar distributions will be automatically resized to ensure this.

Warning

sd_dist will be cloned, rendering it independent of the one passed as input.

compute_corrbool, default=True

If True, returns three values: the Cholesky decomposition, the correlations and the standard deviations of the covariance matrix. Otherwise, only returns the packed Cholesky decomposition. Defaults to True. compatibility.

store_in_tracebool, default=True

Whether to store the correlations and standard deviations of the covariance matrix in the posterior trace. If True, they will automatically be named as {name}_corr and {name}_stds respectively. Effective only when compute_corr=True.

Returns:
cholTensorVariable

If compute_corr=True. The unpacked Cholesky covariance decomposition.

corrTensorVariable

If compute_corr=True. The correlations of the covariance matrix.

stdsTensorVariable

If compute_corr=True. The standard deviations of the covariance matrix.

packed_cholTensorVariable

If compute_corr=False The packed Cholesky covariance decomposition.

Notes

Since the Cholesky factor is a lower triangular matrix, we use packed storage for the matrix: We store the values of the lower triangular matrix in a one-dimensional array, numbered by row:

[[0 - - -]
 [1 2 - -]
 [3 4 5 -]
 [6 7 8 9]]

The unpacked Cholesky covariance matrix is automatically computed and returned when you specify compute_corr=True in pm.LKJCholeskyCov (see example below). Otherwise, you can use pm.expand_packed_triangular(packed_cov, lower=True) to convert the packed Cholesky matrix to a regular two-dimensional array.

References

[1]

Lewandowski, D., Kurowicka, D. and Joe, H. (2009). “Generating random correlation matrices based on vines and extended onion method.” Journal of multivariate analysis, 100(9), pp.1989-2001.

[2]

J. M. isn’t a mathematician (http://math.stackexchange.com/users/498/ j-m-isnt-a-mathematician), Different approaches to evaluate this determinant, URL (version: 2012-04-14): http://math.stackexchange.com/q/130026

Examples

with pm.Model() as model:
    # Note that we access the distribution for the standard
    # deviations, and do not create a new random variable.
    sd_dist = pm.Exponential.dist(1.0, size=10)
    chol, corr, sigmas = pm.LKJCholeskyCov(
        'chol_cov', eta=4, n=10, sd_dist=sd_dist
    )

    # if you only want the packed Cholesky:
    # packed_chol = pm.LKJCholeskyCov(
        'chol_cov', eta=4, n=10, sd_dist=sd_dist, compute_corr=False
    )
    # chol = pm.expand_packed_triangular(10, packed_chol, lower=True)

    # Define a new MvNormal with the given covariance
    vals = pm.MvNormal('vals', mu=np.zeros(10), chol=chol, shape=10)

    # Or transform an uncorrelated normal:
    vals_raw = pm.Normal('vals_raw', mu=0, sigma=1, shape=10)
    vals = pt.dot(chol, vals_raw)

    # Or compute the covariance matrix
    cov = pt.dot(chol, chol.T)

Implementation In the unconstrained space all values of the cholesky factor are stored untransformed, except for the diagonal entries, where we use a log-transform to restrict them to positive values.

To correctly compute log-likelihoods for the standard deviations and the correlation matrix separately, we need to consider a second transformation: Given a cholesky factorization \(LL^T = \Sigma\) of a covariance matrix we can recover the standard deviations \(\sigma\) as the euclidean lengths of the rows of \(L\), and the cholesky factor of the correlation matrix as \(U = \text{diag}(\sigma)^{-1}L\). Since each row of \(U\) has length 1, we do not need to store the diagonal. We define a transformation \(\phi\) such that \(\phi(L)\) is the lower triangular matrix containing the standard deviations \(\sigma\) on the diagonal and the correlation matrix \(U\) below. In this form we can easily compute the different likelihoods separately, as the likelihood of the correlation matrix only depends on the values below the diagonal, and the likelihood of the standard deviation depends only on the diagonal values.

We still need the determinant of the jacobian of \(\phi^{-1}\). If we think of \(\phi\) as an automorphism on \(\mathbb{R}^{\tfrac{n(n+1)}{2}}\), where we order the dimensions as described in the notes above, the jacobian is a block-diagonal matrix, where each block corresponds to one row of \(U\). Each block has arrowhead shape, and we can compute the determinant of that as described in [2]. Since the determinant of a block-diagonal matrix is the product of the determinants of the blocks, we get

\[\text{det}(J_{\phi^{-1}}(U)) = \left[ \prod_{i=2}^N u_{ii}^{i - 1} L_{ii} \right]^{-1}\]

Methods

LKJCholeskyCov.dist(eta, n, sd_dist, *[, ...])