10. Gaussian copulas

The copula submodule provides classes to define probability distributions and reparametrize them such that the joint distribution is Normal. This is useful to define priors for the hyperparameters of a Gaussian process, but can also be used on its own.

To define a variable, use one of the subclasses of Distr listed under Predefined families. Combine the variables together to define a model by putting them in a Copula object. See Distr for examples.

To represent at once concrete values of the variables and their transformed parametrization, put them in a gvar.BufferDict using makedict. To apply the transformation manually, use partial_invfcn.

Note

I define “Gaussian copula” to mean a representation of an arbitrary random variable as the transformation of a multivariate Normal random variable, as explained, e.g., here. This is different from another common usage, which is representing a bivariate Normal as the transformation of uniform variables, as “Gaussian copula” on wikipedia.

10.1. Generic classes

The class DistrBase defines basic functionality shared by Distr and Copula, and can not be instantiated. Distr represents a probability distribution on a numerical tensor; it can not be instantied, use its concrete subclasses listed below. Copula represents a collection of possibly related Distr objects and is intended for direct use.

class lsqfitgp.copula.DistrBase

Abstract base class to represent (trees of) probability distributions.

See also

Distr, Copula

Attributes:

in_shapetuple of int

The core shape of the input array to partial_invfcn.

shape(tree of) tuple of int

The core shape of the output array of partial_invfcn.

dtype(tree of) dtype

The dtype of the output array of partial_invfcn.

distrshape(tree of) tuple of int

The sub-core shape of the output array of partial_invfcn that represents the atomic shape of the distribution.

Methods

partial_invfcn(x)	Map independent Normal variables to the desired distribution.
add_distribution(name)	Register the distribution for usage with gvar.BufferDict.
gvars()	Return an array of gvars intended as value in a gvar.BufferDict.

add_distribution(name)

Register the distribution for usage with gvar.BufferDict.

Parameters:

namestr

The name to use for the distribution. It must be globally unique, and it should not contain parentheses. To redefine a distribution with the same name, use gvar.BufferDict.del_distribution first. However, it is allowed to reuse the name if the distribution family, shape and parameters are identical to those used for the existing definition.

See also

gvar.BufferDict.add_distribution, gvar.BufferDict.del_distribution

gvars()

Return an array of gvars intended as value in a gvar.BufferDict.

Returns:

gvarsarray of gvars

An array of i.i.d. standard Normal primary gvars with shape in_shape.

partial_invfcn(x)

Map independent Normal variables to the desired distribution.

This function is a generalized ufunc. It is jax traceable and differentiable one time. It supports arrays of gvars as input. The attributes in_shape and shape give the core shapes.

Parameters:

x(..., *in_shape) array

An array of values representing draws of i.i.d. Normal variates.

Returns:

y(tree of) (..., *shape) array

An array of values representing draws of the desired distribution.

class lsqfitgp.copula.Distr(*params, shape=(), name=None)

Abstract base class to represent probability distributions.

A Distr object represents a probability distribution of a variable in \(\mathbb R^n\), and provides a transformation function from a (multivariate) Normal variable to the target random variable.

The main functionality is defined in DistrBase. The additional attributes and methods params, signature, and invfcn are not intended for common usage.

Parameters:

*paramstuple of scalar, array or Distr

The parameters of the distribution. If the parameters have leading axes other than those required, the distribution is repeated i.i.d. over those axes. If a parameter is an instance of Distr itself, it is a random parameter and its distribution is accounted for.

shapeint or tuple of int

The shape of the array of i.i.d. variables to be represented, scalar by default. If the variable is multivariate, this shape adds as leading axes in the array. This shape broadcasts with the non-core shapes of the parameters.

namestr, optional

If specified, the distribution is defined for usage with gvar.BufferDict using gvar.BufferDict.add_distribution, and for convenience the constructor returns an array of gvars with the appropriate shape instead of the Distr object. See add_distribution.

Returns:

If name is None (default):

distrDistr

An object representing the distribution.

Else:

gvarsarray of gvars

An array of primary gvars that can be set as value in a gvar.BufferDict under a key that uses the just defined name.

See also

DistrBase, Copula, gvar.BufferDict.uniform

Notes

Concrete subclasses must define invfcn, and define the class attribute signature to the numpy signature string of invfcn, unless invfcn is an ufunc and its number of parameters can be inferred. invfcn must be vectorized.

Examples

Use directly with gvar.BufferDict by setting name:

>>> copula = gvar.BufferDict({
...     'A(x)': lgp.copula.beta(1, 1, name='A'),
...     'B(y)': lgp.copula.beta(3, 5, name='B'),
... })
>>> copula['x']
0.50(40)
>>> copula['y']
0.36(18)

Corresponding “unrolled” usage:

>>> A = lgp.copula.beta(1, 1)
>>> B = lgp.copula.beta(3, 5)
>>> A.add_distribution('A')
>>> B.add_distribution('B')
>>> copula = gvar.BufferDict({
...     'A(x)': A.gvars(),
...     'B(y)': B.gvars(),
... })

Notice that, although the name used for add_distribution must be globally unique, for convenience it is permitted to redefine the same distribution family with the same parameters, even from another Distr instance.

To generate automatically sensible names and avoid repeating them twice, use makedict:

>>> lgp.copula.makedict({
...     'x': lgp.copula.beta(1, 1),
...     'y': lgp.copula.beta(3, 5),
... })
BufferDict({'__copula_beta{1, 1}(x)': 0.0(1.0), '__copula_beta{3, 5}(y)': 0.0(1.0)})

Define a distribution with a random parameter:

>>> X = lgp.copula.halfnorm(np.sqrt(lgp.copula.invgamma(1, 1)))
>>> X
halfnorm(sqrt(invgamma(1, 1)))

Now X represents the model

\[\begin{split}\sigma^2 &\sim \mathrm{InvGamma}(1, 1), \\ X \mid \sigma &\sim \mathrm{HalfNorm}(\sigma).\end{split}\]

In general it is possible to transform a Distr with numpy ufuncs and continuous arithmetic operations.

Repeated usage of Distr instances for random parameters will share those parameters in the distributions. The following code:

>>> sigma2 = lgp.copula.invgamma(1, 1)
>>> X = lgp.copula.halfnorm(np.sqrt(sigma2))
>>> Y = lgp.copula.halfcauchy(np.sqrt(sigma2))

Corresponds to the model

\[\begin{split}\sigma^2 &\sim \mathrm{InvGamma}(1, 1), \\ X \mid \sigma &\sim \mathrm{HalfNorm}(\sigma), \\ Y \mid \sigma &\sim \mathrm{HalfCauchy}(\sigma),\end{split}\]

with the same parameter \(\sigma^2\) shared between the two distributions. However, if the distributions are now put into a gvar.BufferDict, with

>>> sigma2.add_distribution('distr_sigma2')
>>> X.add_distribution('distr_X')
>>> Y.add_distribution('distr_Y')
>>> bd = gvar.BufferDict({
...     'distr_sigma2(sigma2)': sigma2.gvars(),
...     'distr_X(X)': X.gvars(),
...     'distr_Y(Y)': Y.gvars(),
... })

then this relationship breaks down; the model represented by the dictionary bd is

\[\begin{split}\sigma^2 &\sim \mathrm{InvGamma}(1, 1), \\ X \mid \sigma_X &\sim \mathrm{HalfNorm}(\sigma_X), \quad & \sigma_X^2 &\sim \mathrm{InvGamma}(1, 1), \\ Y \mid \sigma_Y &\sim \mathrm{HalfCauchy}(\sigma_Y), \quad & \sigma_Y^2 &\sim \mathrm{InvGamma}(1, 1),\end{split}\]

with separate, independent parameters \(\sigma,\sigma_X,\sigma_Y\), because each dictionary entry is evaluated separately. Indeed, trying to do this with makedict will raise an error:

>>> bd = lgp.copula.makedict({'sigma2': sigma2, 'X': X, 'Y': Y})
ValueError: cross-key occurrences of object(s):
invgamma with id 6201535248: <sigma2>, <X.0.0>, <Y.0.0>

To use all the distributions at once while preserving the relationships, put them into a container of choice and wrap it as a Copula object:

>>> sigmaXY = lgp.copula.Copula({'sigma2': sigma2, 'X': X, 'Y': Y})

The Copula provides a partial_invfcn function to map Normal variables to a structure, with the same layout as the input one, of desired variates. The whole Copula can be used in gvar.BufferDict:

>>> bd = lgp.copula.makedict({'sigmaXY': sigmaXY})
>>> bd
BufferDict({"__copula_{'sigma2': invgamma{1, 1}, 'X': halfnorm{sqrt{_Path{path=[{DictKey{key='sigma2'},}]}}}, 'Y': halfcauchy{sqrt{_Path{path=[{DictKey{key='sigma2'},}]}}}}(sigmaXY)": array([0.0(1.0), 0.0(1.0), 0.0(1.0)], dtype=object)})
>>> bd['sigmaXY']
{'sigma2': 1.4(1.7), 'X': 0.81(89), 'Y': 1.2(1.7)}
>>> gvar.corr(bd['sigmaXY']['X'], bd['sigmaXY']['Y'])
0.21950577757757836

Although the actual dictionary value is a flat array, getting the unwrapped key reproduces the original structure.

To apply arbitrary transformations, use manually invfcn:

>>> @functools.partial(lgp.gvar_gufunc, signature='(n)->(n)')
>>> @functools.partial(jnp.vectorize, signature='(n)->(n)')
>>> def model_invfcn(normal_params):
...     sigma2 = lgp.copula.invgamma.invfcn(normal_params[0], 1, 1)
...     sigma = jnp.sqrt(sigma2)
...     X = lgp.copula.halfnorm.invfcn(normal_params[1], sigma)
...     Y = lgp.copula.halfcauchy.invfcn(normal_params[2], sigma)
...     return jnp.stack([sigma, X, Y])

The jax.numpy.vectorize decorator makes model_invfcn support broadcasting on additional input axes, while gvar_gufunc makes it accept gvars as input.

Attributes:

paramstuple

The parameters as passed to the constructor.

signatureSignature

An object representing the signature of invfcn. This is a class attribute.

Methods

invfcn(x, *params)

Normal to desired distribution transformation.

abstract classmethod invfcn(x, *params)

Normal to desired distribution transformation.

Maps a (multivariate) Normal variable to a variable with the desired marginal distribution. In symbols: \(y = F^{-1}(\Phi(x))\). This function is a generalized ufunc, jax traceable, vmappable one time, and differentiable one time. The signature is accessible through the class attribute signature.

Parameters:

xarray_like

The input Normal variable.

*paramsarray_like

The parameters of the distribution.

Returns:

yarray_like

The output variable with the desired marginal distribution.

class lsqfitgp.copula.Copula(variables)

Represents a tree of probability distributions.

By “tree” it is intended an arbitrarily nested structure of containers (e.g., dict or list) where each leaf node is either a Distr object or another Copula.

The function of a Copula is to keep into account the relationships between all the Distr objects when defining the function partial_invfcn that maps Normal variates to the desired random variable. The same Distr object appearing in multiple places will not be sampled more than once.

This class inherits the functionality defined in DistrBase without additions. The attributes shape, dtype, distrshape, that represents properties of the output of partial_invfcn, are trees matching the structure of the tree of variables. The same holds for the output of partial_invfcn. in_shape instead is an ordinary array shape, indicating that the input Normal variables are a raveled single array.

Parameters:

variablestree of Distr or Copula

The tree of distributions to wrap. The containers are copied, so successive modifications will not be reflected on the Copula.

See also

DistrBase, Distr, jax.tree_util

Examples

Define a model into a dictionary one piece at a time, then wrap it as a Copula:

>>> m = {}
>>> m['a'] = lgp.copula.halfnorm(1.5)
>>> m['b'] = lgp.copula.halfcauchy(2)
>>> m['c'] = [
...     lgp.copula.uniform(m['a'], m['b']),
...     lgp.copula.uniform(m['b'], m['a']),
... ]
>>> cop = lgp.copula.Copula(m)
>>> cop
Copula({'a': halfnorm(1.5),
 'b': halfcauchy(2),
 'c': [uniform(<a>, <b>), uniform(<b>, <a>)]})

Notice how, when showing the object on the REPL, multiple appearances of the same variables are replaced by identifiers derived from the dictionary keys.

The model may then be extended to create a variant. This does not affect cop:

>>> m['d'] = lgp.copula.invgamma(m['c'][0], m['c'][1])
>>> cop2 = lgp.copula.Copula(m)
>>> cop2
Copula({'a': halfnorm(1.5),
 'b': halfcauchy(2),
 'c': [uniform(<a>, <b>), uniform(<b>, <a>)],
 'd': invgamma(<c.0>, <c.1>)})
>>> cop
Copula({'a': halfnorm(1.5),
 'b': halfcauchy(2),
 'c': [uniform(<a>, <b>), uniform(<b>, <a>)]})

10.2. Utilities

lsqfitgp.copula.makedict(variables, prefix='__copula_')

Expand distributions in a dictionary.

Parameters:

variablesdict

A dictionary representing a collection of probability distribution. If a value is an instance of DistrBase, the key is converted to mark a transformation and the value is replaced with new primary gvars.

prefixstr

A prefix to make the transformation names unique.

Returns:

outBufferDict

The transformed dictionary. Recognizes the same keys as variables, but squashes the values through the transformation that sends a Normal to the desired distribution.

Raises:

ValueError

If any DistrBase object appears under different keys.

Examples

Put a Distr into a gvar.BufferDict:

>>> bd = lgp.copula.makedict({'x': lgp.copula.beta(1, 1)})
>>> bd
BufferDict({'__copula_beta{1, 1}(x)': 0.0(1.0)})
>>> bd['x']
0.50(40)
>>> bd['__copula_beta{1, 1}(x)']
0.0(1.0)

You can also put an entire Copula:

>>> bd = lgp.copula.makedict({
...     'x': lgp.copula.Copula({
...         'y': lgp.copula.gamma(2, 1 / 2),
...         'z': lgp.copula.invgamma(2, 2),
...     }),
... })
>>> bd
BufferDict({"__copula_{'y': gamma{2, 0.5}, 'z': invgamma{2, 2}}(x)": array([0.0(1.0), 0.0(1.0)], dtype=object)})
>>> bd['x']
{'y': 3.4(2.5), 'z': 1.19(90)}

Other entries and transformations are left alone:

>>> bd = lgp.copula.makedict({
...     'x': gvar.gvar(3, 0.2),
...     'log(y)': gvar.gvar(0, 1),
...     'z': lgp.copula.dirichlet(1.5, [1, 2, 3]),
... })
>>> bd
BufferDict({'x': 3.00(20), 'log(y)': 0.0(1.0), '__copula_dirichlet{1.5, [1, 2, 3], shape=3}(z)': array([0.0(1.0), 0.0(1.0), 0.0(1.0)], dtype=object)})
>>> bd['z']
array([0.06(20), 0.31(49), 0.63(51)], dtype=object)
>>> bd['y']
1.0(1.0)

Since shared DistrBase objects represent statistical dependency, it is forbidden to have the same object appear under different keys, as that would make it impossible to take the dependencies into account:

>>> x = lgp.copula.beta(1, 1)
>>> y = lgp.copula.beta(1, x)
>>> lgp.copula.makedict({'x': x, 'y': y})
ValueError: cross-key occurrences of object(s):
beta with id 10952248976: <x>, <y.1>

lsqfitgp.copula.distribution(invfcn, signature=None, dtype=None)

Decorator to define a distribution from a transformation function.

Parameters:

invfcnfunction

The transformation function from a (multivariate) standard Normal variable to the target random variable. The signature must be invfcn(x, *params). It must be jax-traceable. It does not need to be vectorized.

signaturestr, optional

The signature of invfcn, as a numpy signature string. If not specified, invfcn is assumed to take and output scalars.

dtypedtype, optional

The dtype of the output of invfcn. If not specified, it is assumed to be floating point.

Returns:

clsDistr

The new distribution class.

Examples

>>> @lgp.copula.distribution
... def uniform(x, a, b):
...     return a + (b - a) * jax.scipy.stats.norm.cdf(x)

>>> @functools.partial(lgp.copula.distribution, signature='(n,m)->(n)')
... def wishart(x):
...     return x @ x.T

10.3. Predefined families

The parametrizations follow Wikipedia, while the class names are as in scipy.stats.

class lsqfitgp.copula.beta(alpha, beta)

https://en.wikipedia.org/wiki/Beta_distribution

class lsqfitgp.copula.dirichlet(alpha, n)

https://en.wikipedia.org/wiki/Dirichlet_distribution

This uses the \((\alpha, \mathbf n)\) parametrization, such that the standard one is

\[\boldsymbol\alpha = \frac{\alpha \mathbf n}{\sum_i n_i}.\]

class lsqfitgp.copula.gamma(alpha, beta)

https://en.wikipedia.org/wiki/Gamma_distribution

class lsqfitgp.copula.halfcauchy(gamma)

https://en.wikipedia.org/wiki/Cauchy_distribution, scipy.stats.halfcauchy

class lsqfitgp.copula.halfnorm(sigma)

https://en.wikipedia.org/wiki/Half-normal_distribution

class lsqfitgp.copula.invgamma(alpha, beta)

https://en.wikipedia.org/wiki/Inverse-gamma_distribution

class lsqfitgp.copula.loggamma(alpha)

https://en.wikipedia.org/wiki/Gamma_distribution, scipy.stats.loggamma

This is the distribution of the logarithm of a Gamma variable. The naming convention is the opposite of lognorm, which is the distribution of the exponential of a Normal variable.

class lsqfitgp.copula.uniform(a, b)

https://en.wikipedia.org/wiki/Continuous_uniform_distribution