gvar - Gaussian Random Variables

Introduction

Objects of type gvar.GVar represent gaussian random variables, which are specified by a mean and standard deviation. They are created using gvar.gvar(): for example,

>>> x = gvar.gvar(10,3)          # 0 +- 3
>>> y = gvar.gvar(12,4)          # 2 +- 4
>>> z = x + y                    # 2 +- 5
>>> print(z)
22.0(5.0)
>>> print(z.mean)
22.0
>>> print(z.sdev)
5.0

This module contains a variety of tools for creating and manipulating gaussian random variables, including:

Functions

The function used to create Gaussian variable objects is:

gvar.gvar(...)

Create one or more new gvar.GVars.

Each of the following creates new gvar.GVars:

gvar.gvar(x, xsdev)

Returns a gvar.GVar with mean x and standard deviation xsdev. Returns an array of gvar.GVars if x and xsdev are arrays with the same shape; the shape of the result is the same as the shape of x. Returns a gvar.BufferDict if x and xsdev are dictionaries with the same keys and layout; the result has the same keys and layout as x.

gvar.gvar(x, xcov)

Returns an array of gvar.GVars with means given by array x and a covariance matrix given by array xcov, where xcov.shape = 2*x.shape; the result has the same shape as x. Returns a gvar.BufferDict if x and xcov are dictionaries, where the keys in xcov are (k1,k2) for any keys k1 and k2 in x. Returns a single gvar.GVar if x is a number and xcov is a one-by-one matrix. The layout for xcov is compatible with that produced by gvar.evalcov() for a single gvar.GVar, an array of gvar.GVars, or a dictionary whose values are gvar.GVars and/or arrays of gvar.GVars. Therefore gvar.gvar(gvar.mean(g), gvar.evalcov(g)) creates gvar.GVars with the same means and covariance matrix as the gvar.GVars in g provided g is a single gvar.GVar, or an array or dictionary of gvar.GVars.

gvar.gvar((x, xsdev))

Returns a gvar.GVar with mean x and standard deviation xsdev.

gvar.gvar(xstr)

Returns a gvar.GVar corresponding to string xstr which is either of the form "xmean +- xsdev" or "x(xerr)" (see GVar.fmt()).

gvar.gvar(xgvar)

Returns gvar.GVar xgvar unchanged.

gvar.gvar(xdict)

Returns a dictionary (BufferDict) b where b[k] = gvar(xdict[k]) for every key in dictionary xdict. The values in xdict, therefore, can be strings, tuples or gvar.GVars (see above), or arrays of these.

gvar.gvar(xarray)

Returns an array a having the same shape as xarray where every element a[i...] = gvar(xarray[i...]). The values in xarray, therefore, can be strings, tuples or gvar.GVars (see above).

gvar.gvar is actually an object of type gvar.GVarFactory.

The following function is useful for constructing new functions that can accept gvar.GVars as arguments:

gvar.gvar_function(x, f, dfdx)

Create a gvar.GVar for function f(x) given f and df/dx at x.

This function creates a gvar.GVar corresponding to a function of gvar.GVars x whose value is f and whose derivatives with respect to each x are given by dfdx. Here x can be a single gvar.GVar, an array of gvar.GVars (for a multidimensional function), or a dictionary whose values are gvar.GVars or arrays of gvar.GVars, while dfdx must be a float, an array of floats, or a dictionary whose values are floats or arrays of floats, respectively.

This function is useful for creating functions that can accept gvar.GVars as arguments. For example,

import math
import gvar as gv

def sin(x):
    if isinstance(x, gv.GVar):
        f = math.sin(x.mean)
        dfdx = math.cos(x.mean)
        return gv.gvar_function(x, f, dfdx)
    else:
        return math.sin(x)

creates a version of sin(x) that works with either floats or gvar.GVars as its argument. This particular function is unnecessary since it is already provided by gvar.

Parameters

  • x (gvar.GVar, array of gvar.GVars, or a dictionary of gvar.GVars) – Point at which the function is evaluated.

  • f (float) – Value of function at point gvar.mean(x).

  • dfdx (float, array of floats, or a dictionary of floats) – Derivatives of function with respect to x at point gvar.mean(x).

Returns

A gvar.GVar representing the function’s value at x.

Means, standard deviations, variances, formatted strings, covariance matrices and correlation/comparison information can be extracted from arrays (or dictionaries) of gvar.GVars using:

gvar.mean(g)

Extract means from gvar.GVars in g.

g can be a gvar.GVar, an array of gvar.GVars, or a dictionary containing gvar.GVars or arrays of gvar.GVars. Result has the same layout as g.

Elements of g that are not gvar.GVars are left unchanged.

gvar.sdev(g)

Extract standard deviations from gvar.GVars in g.

g can be a gvar.GVar, an array of gvar.GVars, or a dictionary containing gvar.GVars or arrays of gvar.GVars. Result has the same layout as g.

The deviation is set to 0.0 for elements of g that are not gvar.GVars.

gvar.var(g)

Extract variances from gvar.GVars in g.

g can be a gvar.GVar, an array of gvar.GVars, or a dictionary containing gvar.GVars or arrays of gvar.GVars. Result has the same layout as g.

The variance is set to 0.0 for elements of g that are not gvar.GVars.

gvar.fmt(g, ndecimal=None, sep='')

Format gvar.GVars in g.

g can be a gvar.GVar, an array of gvar.GVars, or a dictionary containing gvar.GVars or arrays of gvar.GVars. Each gvar.GVar gi in g is replaced by the string generated by gi.fmt(ndecimal,sep). Result has same structure as g.

gvar.tabulate(g, ncol=1, headers=True, offset='', ndecimal=None)

Tabulate contents of an array or dictionary of gvar.GVars.

Given an array g of gvar.GVars or a dictionary whose values are gvar.GVars or arrays of gvar.GVars, gvar.tabulate(g) returns a string containing a table of the values of g’s entries. For example, the code

import collections
import gvar as gv

g = collections.OrderedDict()
g['scalar'] = gv.gvar('10.3(1)')
g['vector'] = gv.gvar(['0.52(3)', '0.09(10)', '1.2(1)'])
g['tensor'] = gv.gvar([
    ['0.01(50)', '0.001(20)', '0.033(15)'],
    ['0.001(20)', '2.00(5)', '0.12(52)'],
    ['0.007(45)', '0.237(4)', '10.23(75)'],
    ])
print(gv.tabulate(g, ncol=2))

prints the following table:

   key/index          value     key/index          value
---------------------------  ---------------------------
      scalar     10.30 (10)           1,0     0.001 (20)
    vector 0     0.520 (30)           1,1     2.000 (50)
           1      0.09 (10)           1,2      0.12 (52)
           2      1.20 (10)           2,0     0.007 (45)
  tensor 0,0      0.01 (50)           2,1    0.2370 (40)
         0,1     0.001 (20)           2,2     10.23 (75)
         0,2     0.033 (15)
Parameters

  • g – Array of gvar.GVars (any shape) or dictionary whose values are gvar.GVars or arrays of gvar.GVars (any shape).

  • ncol – The table is split over ncol columns of key/index values plus gvar.GVar values. Default value is 1.

  • headers – Prints standard header on table if True; omits the header if False. If headers is a 2-tuple, then headers[0] is used in the header over the indices/keys and headers[1] over the gvar.GVar values. (Default is True.)

  • offset (str) – String inserted at the beginning of each line in the table. Default is ''.

  • ndecimal – Number of digits displayed after the decimal point. Default is ndecimal=None which adjusts table entries to show 2 digits of error.

gvar.correlate(g, corr)

Add correlations to uncorrelated gvar.GVars in g.

This method creates correlated gvar.GVars from uncorrelated gvar.GVars g, using the correlations specified in corr.

Note that correlations initially present in g, if any, are ignored.

Examples

A typical application involves the construction of correlated gvar.GVars give the means and standard deviations, together with a correlation matrix:

>>> import gvar as gv
>>> g = gv.gvar(['1(1)', '2(10)'])
>>> print(gv.evalcorr(g))           # uncorrelated
[[ 1.  0.]
 [ 0.  1.]]
>>> g =  gv.correlate(g, [[1. , 0.1], [0.1, 1.]])
>>> print(gv.evalcorr(g))           # correlated
[[ 1.   0.1]
 [ 0.1  1. ]]

This also works when g and corr are dictionaries:

>>> g = gv.gvar(dict(a='1(1)', b='2(10)'))
>>> print(gv.evalcorr(g))
{('a', 'a'): array([[ 1.]]),('a', 'b'): array([[ 0.]]),('b', 'a'): array([[ 0.]]),('b', 'b'): array([[ 1.]])}
>>> corr = {}
>>> corr['a', 'a'] = 1.0
>>> corr['a', 'b'] = 0.1
>>> corr['b', 'a'] = 0.1
>>> corr['b', 'b'] = 1.0
>>> g = correlate(g, corr)
>>> print(gv.evalcorr(g))
{('a', 'a'): array([[ 1.]]),('a', 'b'): array([[ 0.1]]),('b', 'a'): array([[ 0.1]]),('b', 'b'): array([[ 1.]])}
Parameters

  • g – An array of gvar.GVars or a dictionary whose values are gvar.GVars or arrays of gvar.GVars.

  • corr – Correlations between gvar.GVars: corr[i, j] is the correlation between g[i] and g[j].

gvar.evalcov(g)

Compute covariance matrix for elements of array/dictionary g.

If g is an array of gvar.GVars, evalcov returns the covariance matrix as an array with shape g.shape+g.shape. If g is a dictionary whose values are gvar.GVars or arrays of gvar.GVars, the result is a doubly-indexed dictionary where cov[k1,k2] is the covariance for g[k1] and g[k2].

gvar.cov(g1, g2)

Covariance of gvar.GVar g1 with g2.

gvar.evalcov_blocks(g)

Evaluate covariance matrix for elements of g.

Evaluates the covariance matrices for gvar.GVars stored in array or dictionary of arrays/gvar.GVars g. The covariance matrix is decomposed into its block diagonal components, and a list of tuples (idx,bcov) is returned where bcov is a diagonal block of the covariance matrix and idx an array containing the corresponding indices in g.flat for that block. So to reassemble the blocks into a single matrix cov, for example, one would use:

import numpy as np
cov = np.empty((len(g), len(g)), float)
for idx, bcov in evalcov_block(g):
    cov[idx[:, None], idx] = bcov

gvar.evalcov_blocks() is particularly useful when the covariance matrix is sparse; only nonzero elements are retained.

gvar.evalcorr(g)

Compute correlation matrix for elements of array/dictionary g.

If g is an array of gvar.GVars, evalcorr returns the correlation matrix as an array with shape g.shape+g.shape. If g is a dictionary whose values are gvar.GVars or arrays of gvar.GVars, the result is a doubly-indexed dictionary where corr[k1,k2] is the correlation for g[k1] and g[k2].

The correlation matrix is related to the covariance matrix by:

corr[i,j] = cov[i,j] / (cov[i,i] * cov[j,j]) ** 0.5
gvar.corr(g1, g2)

Correlation between gvar.GVars g1 and g2.

gvar.uncorrelated(g1, g2)

Return True if gvar.GVars in g1 uncorrelated with those in g2.

g1 and g2 can be gvar.GVars, arrays of gvar.GVars, or dictionaries containing gvar.GVars or arrays of gvar.GVars. Returns True if either of g1 or g2 is None.

gvar.chi2(g1, g2, svdcut=1e-12)

Compute chi**2 of g1-g2.

chi**2 is a measure of whether the multi-dimensional Gaussian distributions g1 and g2 (dictionaries or arrays) agree with each other — that is, do their means agree within errors for corresponding elements. The probability is high if chi2(g1,g2)/chi2.dof is of order 1 or smaller.

Usually g1 and g2 are dictionaries with the same keys, where g1[k] and g2[k] are gvar.GVars or arrays of gvar.GVars having the same shape. Alternatively g1 and g2 can be gvar.GVars, or arrays of gvar.GVars having the same shape.

One of g1 or g2 can contain numbers instead of gvar.GVars, in which case chi**2 is a measure of the likelihood that the numbers came from the distribution specified by the other argument.

One or the other of g1 or g2 can be missing keys, or missing elements from arrays. Only the parts of g1 and g2 that overlap are used. Also setting g2=None is equivalent to replacing its elements by zeros.

chi**2 is computed from the inverse of the covariance matrix of g1-g2. The matrix inversion can be sensitive to roundoff errors. In such cases, SVD cuts can be applied by setting parameters svdcut; see the documentation for gvar.svd(), which is used to apply the cut.

The return value is the chi**2. Extra attributes attached to this value give additional information:

  • dof — Number of degrees of freedom (that is, the number of variables compared).

  • Q — The probability that the chi**2 could have been larger, by chance, even if g1 and g2 agree. Values smaller than 0.1 or so suggest that they do not agree. Also called the p-value.

gvar.fmt_chi2(f)

Return string containing chi**2/dof, dof and Q from f.

Assumes f has attributes chi2, dof and Q. The logarithm of the Bayes factor will also be printed if f has attribute logGBF.

gvar.GVars are compared by:

gvar.equivalent(g1, g2, rtol=1e-10, atol=1e-10)

Determine whether g1 and g2 contain equivalent gvar.GVars.

Compares sums and differences of gvar.GVars stored in g1 and g2 to see if they agree within tolerances. Operationally, agreement means that:

abs(diff) < abs(summ) / 2 * rtol + atol

where diff and summ are the difference and sum of the mean values (g.mean) or derivatives (g.der) associated with each pair of gvar.GVars.

gvar.GVars that are equivalent are effectively interchangeable with respect to both their means and also their covariances with any other gvar.GVar (including ones not in g1 and g2).

g1 and g2 can be individual gvar.GVars or arrays of gvar.GVars or dictionaries whose values are gvar.GVars and/or arrays of gvar.GVars. Comparisons are made only for shared keys when they are dictionaries. Array dimensions must match between g1 and g2, but the shapes can be different; comparisons are made for the parts of the arrays that overlap in shape.

Parameters

  • g1 – A gvar.GVar or an array of gvar.GVars or a dictionary of gvar.GVars and/or arrays of gvar.GVars.

  • g2 – A gvar.GVar or an array of gvar.GVars or a dictionary of gvar.GVars and/or arrays of gvar.GVars.

  • rtol – Relative tolerance with which mean values and derivatives must agree with each other. Default is 1e-10.

  • atol – Absolute tolerance within which mean values and derivatives must agree with each other. Default is 1e-10.

gvar.GVars can be stored (serialized) and retrieved from files (or strings) using:

gvar.dump(g, outputfile, method='pickle')

Serialize a collection g of gvar.GVars into file outputfile.

The gvar.GVars are recovered using gvar.load().

Three serialization methods are available: pickle, json, and yaml (provided the yaml module is installed).

json can have trouble with dictionaries whose keys are not strings. A workaround is used here that succeeds provided eval(repr(k)) == k for every key k, which is true for strings and lots of other types of key. Use pickle where the workaround fails.

Parameters

  • g – A gvar.GVar, array of gvar.GVars, or dictionary whose values are gvar.GVars and/or arrays of gvar.GVars.

  • outputfile – The name of a file or a file object in which the serialized gvar.GVars are stored.

  • method (str) – Serialization method, which should be one of ['pickle', 'json', 'yaml']. Default is 'pickle'.

gvar.dumps(g, method='pickle')

Serialize a collection g of gvar.GVars into a string.

The gvar.GVars are recovered using gvar.loads().

Three serialization methods are available: pickle, json, and yaml (provided the yaml module is installed).

json can have trouble with dictionaries whose keys are not strings. A workaround is used here that succeeds provided eval(repr(k)) == k for every key k, which is true for strings and lots of other types of key. Use pickle where the workaround fails.

Parameters

  • g – A gvar.GVar, array of gvar.GVars, or dictionary whose values are gvar.GVars and/or arrays of gvar.GVars.

  • method (str) – Serialization method, which should be one of ['pickle', 'json', 'yaml']. Default is 'pickle'.

gvar.load(inputfile, method=None)

Load and return serialized gvar.GVars from file inputfile.

This function recovers gvar.GVars pickled with gvar.dump().

Parameters

  • inputfile – The name of the file or a file object in which the serialized gvar.GVars are stored.

  • method (str or None) – Serialization method, which should be one of ['pickle', 'json', 'yaml']. If method=None, then each method is tried in turn.

Returns

The reconstructed gvar.GVar, or array or dictionary of gvar.GVars.

gvar.loads(inputstring, method=None)

Load and return serialized gvar.GVars from string inputstring.

This function recovers gvar.GVars pickled with gvar.dumps().

Parameters

  • inputstring – A string containing gvar.GVars serialized using gvar.dumps().

  • method (str or None) – Serialization method, which should be one of ['pickle', 'json', 'yaml']. If method=None, then each method is tried in turn.

Returns

The reconstructed gvar.GVar, or array or dictionary of gvar.GVars.

gvar.disassemble(g)

Disassemble collection g of gvar.GVars.

Disassembles collection g of gvar.GVars into components that can be pickled or otherwise stored. The output is reassembled by gvar.reassemble().

Parameters

g (dict, array, or gvar.GVar) – Collection of gvar.GVars to be disassembled.

gvar.reassemble(data, cov=gvar.gvar.cov)

Convert data (from disassemble) back into gvar.GVars.

Parameters

gvar.GVars contain information about derivatives with respect to the independent gvar.GVars from which they were constructed. This information can be extracted using:

gvar.deriv(g, x)

Compute first derivatives wrt x of gvar.GVars in g.

g can be a gvar.GVar, an array of gvar.GVars, or a dictionary containing gvar.GVars or arrays of gvar.GVars. Result has the same layout as g.

x must be an primary gvar.GVar, which is a gvar.GVar created by a call to gvar.gvar() (e.g., x = gvar.gvar(xmean, xsdev)) or a function f(x) of such a gvar.GVar. (More precisely, x.der must have only one nonzero entry.)

The following functions are used to generate random arrays or dictionaries from the distribution defined by array (or dictionary) g of gvar.GVars. The random numbers incorporate any correlations implied by the gs.

gvar.raniter(g, n=None, svdcut=1e-12)

Return iterator for random samples from distribution g

The gaussian variables (gvar.GVar objects) in array (or dictionary) g collectively define a multidimensional gaussian distribution. The iterator defined by raniter() generates an array (or dictionary) containing random numbers drawn from that distribution, with correlations intact.

The layout for the result is the same as for g. So an array of the same shape is returned if g is an array. When g is a dictionary, individual entries g[k] may be gvar.GVars or arrays of gvar.GVars, with arbitrary shapes.

raniter() also works when g is a single gvar.GVar, in which case the resulting iterator returns random numbers drawn from the distribution specified by g.

Parameters

  • g – An array (or dictionary) of objects of type gvar.GVar; or a gvar.GVar.

  • n (int or None) – Maximum number of random iterations. Setting n=None (the default) implies there is no maximum number.

  • svdcut (float or None) – If positive, replace eigenvalues eig of g’s correlation matrix with max(eig, svdcut * max_eig) where max_eig is the largest eigenvalue; if negative, discard eigenmodes with eigenvalues smaller than |svdcut| * max_eig. Default is 1e-12.

Returns

An iterator that returns random arrays or dictionaries with the same shape as g drawn from the Gaussian distribution defined by g.

gvar.sample(g, svdcut=1e-12)

Generate random sample from distribution g.

Equivalent to next(gvar.raniter(g, svdcut=svdcut)).

Parameters

  • g – An array or dictionary of objects of type gvar.GVar; or a gvar.GVar.

  • svdcut (float or None) – If positive, replace eigenvalues eig of g’s correlation matrix with max(eig, svdcut * max_eig) where max_eig is the largest eigenvalue; if negative, discard eigenmodes with eigenvalues smaller than |svdcut| * max_eig. Default is 1e-12.

Returns

A random array or dictionary, with the same shape as g, drawn from the Gaussian distribution defined by g.

gvar.bootstrap_iter(g, n=None, svdcut=1e-12)

Return iterator for bootstrap copies of g.

The gaussian variables (gvar.GVar objects) in array (or dictionary) g collectively define a multidimensional gaussian distribution. The iterator created by bootstrap_iter() generates an array (or dictionary) of new gvar.GVars whose covariance matrix is the same as g’s but whose means are drawn at random from the original g distribution. This is a bootstrap copy of the original distribution. Each iteration of the iterator has different means (but the same covariance matrix).

bootstrap_iter() also works when g is a single gvar.GVar, in which case the resulting iterator returns bootstrap copies of the g.

Parameters

  • g (array or dictionary or BufferDict) – An array (or dictionary) of objects of type gvar.GVar.

  • n – Maximum number of random iterations. Setting n=None (the default) implies there is no maximum number.

  • svdcut (None or number) – If positive, replace eigenvalues eig of g’s correlation matrix with max(eig, svdcut * max_eig) where max_eig is the largest eigenvalue; if negative, discard eigenmodes with eigenvalues smaller than |svdcut| * max_eig. Default is 1e-12.

Returns

An iterator that returns bootstrap copies of g.

This function is used to seed the random number generator used by gvar:

gvar.ranseed(a)

Seed random number generators with tuple seed.

Argument seed is an integer or a tuple of integers that is used to seed the random number generators used by numpy and random (and therefore by gvar). Reusing the same seed results in the same set of random numbers.

ranseed generates its own seed when called without an argument or with seed=None. This seed is stored in ranseed.seed and also returned by the function. The seed can be used to regenerate the same set of random numbers at a later time.

Parameters

seed (int, tuple, or None) – Seed for generator. Generates a random tuple if None.

Returns

The seed used to reseed the generator.

The following two functions that are useful for tabulating results and for analyzing where the errors in a gvar.GVar constructed from other gvar.GVars come from:

gvar.fmt_errorbudget(outputs, inputs, ndecimal=2, percent=True, verify=False, colwidth=10)

Tabulate error budget for outputs[ko] due to inputs[ki].

For each output outputs[ko], fmt_errorbudget computes the contributions to outputs[ko]’s standard deviation coming from the gvar.GVars collected in inputs[ki]. This is done for each key combination (ko,ki) and the results are tabulated with columns and rows labeled by ko and ki, respectively. If a gvar.GVar in inputs[ki] is correlated with other gvar.GVars, the contribution from the others is included in the ki contribution as well (since contributions from correlated gvar.GVars cannot be distinguished). The table is returned as a string.

Parameters

  • outputs – Dictionary of gvar.GVars for which an error budget is computed.

  • inputs – Dictionary of: gvar.GVars, arrays/dictionaries of gvar.GVars, or lists of gvar.GVars and/or arrays/dictionaries of gvar.GVars. fmt_errorbudget tabulates the parts of the standard deviations of each outputs[ko] due to each inputs[ki].

  • ndecimal (int) – Number of decimal places displayed in table.

  • percent (boolean) – Tabulate % errors if percent is True; otherwise tabulate the errors themselves.

  • colwidth (positive integer or None) – Width of each column. This is set automatically, to accommodate label widths, if colwidth=None (default).

  • verify (boolean) – If True, a warning is issued if: 1) different inputs are correlated (and therefore double count errors); or 2) the sum (in quadrature) of partial errors is not equal to the total error to within 0.1% of the error (and the error budget is incomplete or overcomplete). No checking is done if verify==False (default).

Returns

A table (str) containing the error budget. Output variables are labeled by the keys in outputs (columns); sources of uncertainty are labeled by the keys in inputs (rows).

gvar.fmt_values(outputs, ndecimal=None)

Tabulate gvar.GVars in outputs.

Parameters

  • outputs – A dictionary of gvar.GVar objects.

  • ndecimal (int or None) – Format values v using v.fmt(ndecimal).

Returns

A table (str) containing values and standard deviations for variables in outputs, labeled by the keys in outputs.

The following function applies an SVD cut to the correlation matrix of a set of gvar.GVars:

gvar.svd(g, svdcut=1e-12, wgts=False, add_svdnoise=False)

Apply SVD cuts to collection of gvar.GVars in g.

Standard usage is, for example,

svdcut = ...
gmod = svd(g, svdcut=svdcut)

where g is an array of gvar.GVars or a dictionary containing gvar.GVars and/or arrays of gvar.GVars. When svdcut>0, gmod is a copy of g whose gvar.GVars have been modified to make their correlation matrix less singular than that of the original g: each eigenvalue eig of the correlation matrix is replaced by max(eig, svdcut * max_eig) where max_eig is the largest eigenvalue. This SVD cut, which is applied separately to each block-diagonal sub-matrix of the correlation matrix, increases the variance of the eigenmodes with eigenvalues smaller than svdcut * max_eig.

The modification of g’s covariance matrix is implemented by adding (to g) a set of gvar.GVars with zero means:

gmod = g + gmod.svdcorrection

where gmod.svdcorrection is an array/dictionary containing the gvar.GVars. If parameter add_svdnoise=True, noise is included in gmod.svdcorrection,

gmod.svdcorrection += gv.sample(gmod.svdcorrection),

before it is added to g. The noise can be useful for testing fits and other applications.

When svdcut is negative, eigenmodes of the correlation matrix whose eigenvalues are smaller than |svdcut| * max_eig are dropped from the new matrix and the corresponding components of g are zeroed out (that is, replaced by 0(0)) in gmod.

There is an additional parameter wgts in gvar.svd() whose default value is False. Setting wgts=1 or wgts=-1 instead causes gvar.svd() to return a tuple (gmod, i_wgts) where gmod is the modified copy of g, and i_wgts contains a spectral decomposition of the covariance matrix corresponding to the modified correlation matrix if wgts=1, or a decomposition of its inverse if wgts=-1. The first entry i, wgts = i_wgts[0] specifies the diagonal part of the matrix: i is a list of the indices in gmod.flat corresponding to diagonal elements, and wgts ** 2 gives the corresponding matrix elements. The second and subsequent entries, i, wgts = i_wgts[n] for n > 0, each correspond to block-diagonal sub-matrices, where i is the list of indices corresponding to the block, and wgts[j] are eigenvectors of the sub-matrix rescaled so that

numpy.sum(numpy.outer(wi, wi) for wi in wgts[j]

is the sub-matrix (wgts=1) or its inverse (wgts=-1).

To compute the inverse of the covariance matrix from i_wgts, for example, one could use code like:

gmod, i_wgts = svd(g, svdcut=svdcut, wgts=-1)

inv_cov = numpy.zeros((n, n), float)
i, wgts = i_wgts[0]                       # 1x1 sub-matrices
if len(i) > 0:
    inv_cov[i, i] = numpy.array(wgts) ** 2
for i, wgts in i_wgts[1:]:                # nxn sub-matrices (n>1)
    for w in wgts:
        inv_cov[i[:, None], i] += numpy.outer(w, w)

This sets inv_cov equal to the inverse of the covariance matrix of the gmods. Similarly, we can compute the expectation value, u.dot(inv_cov.dot(v)), between two vectors (numpy arrays) using:

result = 0.0
i, wgts = i_wgts[0]                       # 1x1 sub-matrices
if len(i) > 0:
    result += numpy.sum((u[i] * wgts) * (v[i] * wgts))
for i, wgts in i_wgts[1:]:                # nxn sub-matrices (n>1)
    result += numpy.sum(wgts.dot(u[i]) * wgts.dot(v[i]))

where result is the desired expectation value.

Parameters

  • g – An array of gvar.GVars or a dicitionary whose values are gvar.GVars and/or arrays of gvar.GVars.

  • svdcut (None or float) – If positive, replace eigenvalues eig of the correlation matrix with max(eig, svdcut * max_eig) where max_eig is the largest eigenvalue; if negative, discard eigenmodes with eigenvalues smaller than |svdcut| * max_eig. Note |svdcut| < 1. Default is 1e-12.

  • wgts – Setting wgts=1 causes gvar.svd() to compute and return a spectral decomposition of the covariance matrix of the modified gvar.GVars, gmod. Setting wgts=-1 results in a decomposition of the inverse of the covariance matrix. The default value is False, in which case only gmod is returned.

  • add_svdnoise – If True, noise is added to the SVD correction (see above).

Returns

A copy gmod of g whose correlation matrix is modified by SVD cuts. If wgts is not False, a tuple (g, i_wgts) is returned where i_wgts contains a spectral decomposition of gmod’s covariance matrix or its inverse.

Data from the SVD analysis is stored in gmod:

gmod.svdcut

SVD cut used to create gmod.

gmod.dof

Number of independent degrees of freedom left after the SVD cut. This is the same as the number initially unless svdcut < 0 in which case it may be smaller.

gmod.nmod

Number of modes whose eignevalue was modified by the SVD cut.

gmod.nblocks

A dictionary where gmod.nblocks[s] contains the number of block-diagonal s-by-s sub-matrices in the correlation matrix.

gmod.eigen_range

Ratio of the smallest to largest eigenvalue before SVD cuts are applied (but after rescaling).

gmod.logdet

Logarithm of the determinant of the covariance matrix after SVD cuts are applied (excluding any omitted modes when svdcut < 0 and any diagonal zero modes).

gmod.svdcorrection

Array or dictionary containing the SVD corrections added to g to create gmod: gmod = g + gmod.svdcorrection.

This function is useful when the correlation matrix is singular or almost singular, and its inverse is needed (as in curve fitting).

The following function can be used to rebuild collections of gvar.GVars, ignoring all correlations with other variables. It can also be used to introduce correlations between uncorrelated variables.

gvar.rebuild(g, gvar=gvar, corr=0.0)

Rebuild g stripping correlations with variables not in g.

g is either an array of gvar.GVars or a dictionary containing gvar.GVars and/or arrays of gvar.GVars. rebuild(g) creates a new collection gvar.GVars with the same layout, means and covariance matrix as those in g, but discarding all correlations with variables not in g.

If corr is nonzero, rebuild will introduce correlations wherever there aren’t any using

cov[i,j] -> corr * sqrt(cov[i,i]*cov[j,j])

wherever cov[i,j]==0.0 initially. Positive values for corr introduce positive correlations, negative values anti-correlations.

Parameter gvar specifies a function for creating new gvar.GVars that replaces gvar.gvar() (the default).

Parameters

  • g (array or dictionary) – gvar.GVars to be rebuilt.

  • gvar (gvar.GVarFactory or None) – Replacement for gvar.gvar() to use in rebuilding. Default is gvar.gvar().

  • corr (number) – Size of correlations to introduce where none exist initially.

Returns

Array or dictionary (gvar.BufferDict) of gvar.GVars (same layout as g) where all correlations with variables other than those in g are erased.

The following functions creates new functions that generate gvar.GVars (to replace gvar.gvar()):

gvar.switch_gvar()

Switch gvar.gvar() to new gvar.GVarFactory.

Returns

New gvar.gvar().

gvar.restore_gvar()

Restore previous gvar.gvar().

Returns

Previous gvar.gvar().

gvar.gvar_factory(cov=None)

Return new function for creating gvar.GVars (to replace gvar.gvar()).

If cov is specified, it is used as the covariance matrix for new gvar.GVars created by the function returned by gvar_factory(cov). Otherwise a new covariance matrix is created internally.

gvar.GVars created by different functions cannot be combined in arithmetic expressions (the error message “Incompatible GVars.” results).

gvar.GVar Objects

The fundamental class for representing Gaussian variables is:

class gvar.GVar

The basic attributes are:

mean

Mean value.

sdev

Standard deviation.

var

Variance.

Two methods allow one to isolate the contributions to the variance or standard deviation coming from other gvar.GVars:

partialvar(*args)

Compute partial variance due to gvar.GVars in args.

This method computes the part of self.var due to the gvar.GVars in args. If args[i] is correlated with other gvar.GVars, the variance coming from these is included in the result as well. (This last convention is necessary because variances associated with correlated gvar.GVars cannot be disentangled into contributions corresponding to each variable separately.)

Parameters

args[i] – A gvar.GVar or array/dictionary of gvar.GVars contributing to the partial variance.

Returns

Partial variance due to all of args.

partialsdev(*args)

Compute partial standard deviation due to gvar.GVars in args.

This method computes the part of self.sdev due to the gvar.GVars in args. If args[i] is correlated with other gvar.GVars, the standard deviation coming from these is included in the result as well. (This last convention is necessary because variances associated with correlated gvar.GVars cannot be disentangled into contributions corresponding to each variable separately.)

Parameters

args[i] (gvar.GVar or array/dictionary of gvar.GVars) – Variables contributing to the partial standard deviation.

Returns

Partial standard deviation due to args.

Partial derivatives of the gvar.GVar with respect to the independent gvar.GVars from which it was constructed are given by:

deriv(x)

Derivative of self with respest to primary gvar.GVar x.

All gvar.GVars are constructed from primary gvar.GVars. self.deriv(x) returns the partial derivative of self with respect to primary gvar.GVar x, holding all of the other primary gvar.GVars constant.

Parameters

x – A primary gvar.GVar (or a function of a single primary gvar.GVar).

Returns

The derivative of self with respect to x.

There are two methods for converting self into a string, for printing:

__str__()

Return string representation of self.

The representation is designed to show at least one digit of the mean and two digits of the standard deviation. For cases where mean and standard deviation are not too different in magnitude, the representation is of the form 'mean(sdev)'. When this is not possible, the string has the form 'mean +- sdev'.

fmt(ndecimal=None, sep='')

Convert to string with format: mean(sdev).

Leading zeros in the standard deviation are omitted: for example, 25.67 +- 0.02 becomes 25.67(2). Parameter ndecimal specifies how many digits follow the decimal point in the mean. Parameter sep is a string that is inserted between the mean and the (sdev). If ndecimal is None (default), it is set automatically to the larger of int(2-log10(self.sdev)) or 0; this will display at least two digits of error. Very large or very small numbers are written with exponential notation when ndecimal is None.

Setting ndecimal < 0 returns mean +- sdev.

Two attributes and a method make reference to the original variables from which self is derived:

cov

Underlying covariance matrix (type gvar.smat) shared by all gvar.GVars.

der

Array of derivatives with respect to underlying (original) gvar.GVars.

dotder(v)

Return the dot product of self.der and v.

gvar.BufferDict Objects

gvar.BufferDict objects are ordered dictionaries that are heavily used in gvar’s implementation. They provide the most flexible representation for multi-dimensional Gaussian distributions.

gvar.BufferDict objects differ from ordinary dictionaries in two respects. The first difference is that the dictionary’s values must be scalars or numpy arrays (any shape) of scalars. The scalars can be ordinary integers or floats, but the dictionary was designed especially for gvar.GVars. The dictionary’s values are packed into different parts of a single one-dimensional array or buffer. Items can be added one at a time as in other dictionaries,

>>> import gvar as gv
>>> b = gv.BufferDict()
>>> b['s'] = 0.0
>>> b['v'] = [1., 2.]
>>> b['t'] = [[3., 4.], [5., 6.]]
>>> print(b)
{'s': 0.0,'v': array([1., 2.]),'t': array([[3., 4.],
       [5., 6.]])}

but the values can also be accessed all at once through the buffer:

>>> print(b.buf)
[0. 1. 2. 3. 4. 5. 6.]

The size, shape, and type of the data associated with a given key cannot be changed by later assignments. For example, here b['s'] = [10., 20.] would generate a ValueError exception. The actual value can, of course, be modified: b['s'] = 22. is fine.

The second difference between gvar.BufferDicts and other dictionaries is illustrated by the following code:

>>> b = gv.BufferDict()
>>> b['log(a)'] = gv.gvar('1(1)')
>>> print(b)
{'log(a)': 1.0(1.0)}
>>> print(b['a'], b['log(a)'])
2.7(2.7) 1.0(1.0)

Even though 'a' is not a key in the dictionary, b['a'] is still defined: it equals exp(b['log(a)']). This feature is used to provide (limited) support for non-Gaussian distributions. Here gvar.BufferDict b specifies a distribution that is Gaussain for p['log(a)'], and therefore log-normal for p['a']. Thus, for example,

>>> [x['a'] for x in gv.raniter(b, n=4)]
[2.1662650927997817, 2.3350022125310317, 8.732161128765775, 3.578188553455522]

creates a list of four random numbers drawn from a log-normal distribution. Note that 'a' in this example is not a key in dictionary b, even though both b['a'] and b.get('a') return values:

>>> print('a' in b)
False
>>> print('log(a)' in b)
True
>>> print(list(b))
['log(a)']
>>> print(b['a'], b.get('a'))
2.7(2.7) 2.7(2.7)

This functionality is used routinely by other modules (e.g., lsqfit). The supported distributions, and methods for adding new ones, are described in the documentation for gvar.BufferDict.add_distribution(), below.

class gvar.BufferDict

Ordered dictionary whose data are packed into a 1-d buffer.

gvar.BufferDicts can be created in the usual way dictionaries are created:

>>> b = BufferDict()
>>> b['a'] = 1.
>>> b['b'] = 2
>>> print(b)
{'a': 1.0,'b':2.0}
>>> b = BufferDict(a=1., b=2.)
>>> b = BufferDict([('a',1.), ('b',2.)])

They can also be created from other dictionaries or gvar.BufferDicts:

>>> c = BufferDict(b)
>>> print(c)
{'a': 1.0,'b': 2.0}
>>> c = BufferDict(b, keys=['b'])
>>> print(c)
{'b': 2.0}

The values in a gvar.BufferDict are scalars or arrays of a scalar type (gvar.GVar, float, int, etc.). The data type is normally inferred (dynamically) from the data itself, but can be specified when creating the gvar.BufferDict from another dictionary or list, using keyword dtype:

>>> b = BufferDict(dict(a=1.2), dtype=int)
>>> print(b, b.dtype)
{'a': 1} int64

Some simple arithemetic is allowed between two gvar.BufferDicts, say, g1 and g2 provided they have the same keys and array shapes. So, for example:

>>> a = BufferDict(a=1., b=[2., 3.])
>>> b = BufferDict(a=10., b=[20., 30.])
>>> print(a + b)
    {'a': 11.0,'b': array([22., 33.])}

Subtraction is also defined, as are multiplication and division by scalars. The corresponding +=, -=, *=, /= operators are supported, as are unary + and -.

Finally a gvar.BufferDict can be cloned from another one but with a different buffer (containing different values):

>>> b = BufferDict(a=1., b=2.)
>>> c = BufferDict(b, buf=[10, 20])
>>> print(c)
{'a': 10,'b': 20}

The main attributes are:

size

Size of buffer array.

flat

Buffer array iterator.

dtype

Data type of buffer array elements.

buf

The (1d) buffer array. Allows direct access to the buffer: for example, self.buf[i] = new_val sets the value of the i-th element in the buffer to value new_val. Setting self.buf = nbuf replaces the old buffer by new buffer nbuf. This only works if nbuf is a one-dimensional numpy array having the same length as the old buffer, since nbuf itself is used as the new buffer (not a copy).

shape

Always equal to None. This attribute is included since gvar.BufferDicts share several attributes with numpy arrays to simplify coding that might support either type. Being dictionaries they do not have shapes in the sense of numpy arrays (hence the shape is None).

In addition to standard dictionary methods, the main methods here are:

flatten()

Copy of buffer array.

slice(k)

Return slice/index in self.flat corresponding to key k.

slice_shape(k)

Return tuple (slice/index, shape) corresponding to key k.

has_dictkey(k)

Returns True if self[k] is defined; False otherwise.

Note that k may be a key or it may be related to a related key associated with a non-Gaussian distribution (e.g., 'log(k)'; see gvar.BufferDict.add_distribution() for more information).

static add_distribution(name, invfcn)

Add new parameter distribution.

gvar.BufferDicts can be used to represent a restricted class of non-Gaussian distributions. For example, the code

import gvar as gv
gv.BufferDict.add_distribution('log', gv.exp)

enables the use of log-normal distributions for parameters. So defining, for example,

b = gv.BufferDict()
b['log(a)'] = gv.gvar('1(1)')

means that b['a'] has a value (equal to exp(b['log(a)']) even though 'a' is not a key in the dictionary.

The distributions available by default correspond to:

gv.BufferDict.add_distribution('log', gv.exp)
gv.BufferDict.add_distribution('sqrt', gv.square)
gv.BufferDict.add_distribution('erfinv', gv.erf)

Additional distributions can be added by specifying:

Parameters

  • name (str) – Distributions’ function name.

  • invfcn (callable) – Inverse of the transformation function.

static del_distribution(invfcn)

Delete gvar.BufferDict distribution name.

gvar.SVD Objects

SVD analysis is handled by the following class:

class gvar.SVD(mat, svdcut=None, svdnum=None, compute_delta=False, rescale=False)

SVD decomposition of a pos. sym. matrix.

SVD is a function-class that computes the eigenvalues and eigenvectors of a positive symmetric matrix mat. Eigenvalues that are small (or negative, because of roundoff) can be eliminated or modified using svd cuts. Typical usage is:

>>> mat = [[1.,.25],[.25,2.]]
>>> s = SVD(mat)
>>> print(s.val)             # eigenvalues
[ 0.94098301  2.05901699]
>>> print(s.vec[0])          # 1st eigenvector (for s.val[0])
[ 0.97324899 -0.22975292]
>>> print(s.vec[1])          # 2nd eigenvector (for s.val[1])
[ 0.22975292  0.97324899]

>>> s = SVD(mat,svdcut=0.6)  # force s.val[i]>=s.val[-1]*0.6
>>> print(s.val)
[ 1.2354102   2.05901699]
>>> print(s.vec[0])          # eigenvector unchanged
[ 0.97324899 -0.22975292]

>>> s = SVD(mat)
>>> w = s.decomp(-1)         # decomposition of inverse of mat
>>> invmat = sum(numpy.outer(wj,wj) for wj in w)
>>> print(numpy.dot(mat,invmat))    # should be unit matrix
[[  1.00000000e+00   2.77555756e-17]
 [  1.66533454e-16   1.00000000e+00]]

Input parameters are:

Parameters

  • mat (2-d sequence (numpy.array or list or …)) – Positive, symmetric matrix.

  • svdcut (None or number (|svdcut|<=1).) – If positive, replace eigenvalues of mat with svdcut*(max eigenvalue); if negative, discard eigenmodes with eigenvalues smaller than svdcut times the maximum eigenvalue.

  • svdnum (None or int) – If positive, keep only the modes with the largest svdnum eigenvalues; ignore if set to None.

  • compute_delta (boolean) – Compute delta (see below) if True; set delta=None otherwise.

  • rescale – Rescale the input matrix to make its diagonal elements equal to +-1.0 before diagonalizing.

The results are accessed using:

val

An ordered array containing the eigenvalues of mat. Note that val[i]<=val[i+1].

vec

Eigenvectors vec[i] corresponding to the eigenvalues val[i].

valmin

Minimum eigenvalue allowed in the modified matrix.

valorig

Eigenvalues of original matrix.

D

The diagonal matrix used to precondition the input matrix if rescale==True. The matrix diagonalized is D M D where M is the input matrix. D is stored as a one-dimensional vector of diagonal elements. D is None if rescale==False.

nmod

The first nmod eigenvalues in self.val were modified by the SVD cut (equals 0 unless svdcut > 0).

eigen_range

Ratio of the smallest to the largest eigenvector in the unconditioned matrix (after rescaling if rescale=True)

delta

A vector of gvars whose means are zero and whose covariance matrix is what was added to mat to condition its eigenvalues. Is None if svdcut<0 or compute_delta==False.

decomp(n)

Vector decomposition of input matrix raised to power n.

Computes vectors w[i] such that

mat**n = sum_i numpy.outer(w[i],w[i])

where mat is the original input matrix to svd. This decomposition cannot be computed if the input matrix was rescaled (rescale=True) except for n=1 and n=-1.

Parameters

n (number) – Power of input matrix.

Returns

Array w of vectors.

Requirements

gvar makes heavy use of numpy for array manipulations. It also uses the numpy code for implementing elementary functions (e.g., sin, exp …) in terms of member functions.