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:

• mean(g) — extract means.
• sdev(g) — extract standard deviations.
• var(g) — extract variances.
• fmt(g) — replace all gvar.GVars in array/dictionary by string representations.
• tabulate(g) — tabulate entries in array/dictionary of gvar.GVars.
• correlate(g, corr) — add correlations to gvar.GVars in array/dictionary g.
• chi2(g1, g2) — chi**2 of g1-g2.
• equivalent(g1, g2) — gvar.GVars the same in g1 and g2?
• evalcov(g) — compute covariance matrix.
• evalcov_blocks(g) — compute diagonal blocks of covariance matrix.
• evalcorr(g) — compute correlation matrix.
• fmt_values(g) — list values for printing.
• fmt_errorbudget(g) — create error-budget table for printing.
• fmt_chi2(f) — format chi**2 information in f as string for printing.
• class BufferDict — ordered dictionary with data buffer.
• class PDF — probability density function.
• class PDFStatistics — statistical analysis of moments of a random variable.
• class PDFHistogram — tool for building PDF histograms.
• dump(g, outputfile) — serialize a collection of gvar.GVars in file.
• dumps(g) — serialize a collection of gvar.GVars in a string.
• load(inputfile) — reconstitute a collection of gvar.GVars from a file.
• loads(inputstr) — reconstitute a collection of gvar.GVars from a string.
• disassemble(g) — low-level routine to disassemble a collection of gvar.GVars.
• reassemble(data,cov) — low-level routine to reassemble a collection of gvar.GVars.
• raniter(g,N) — iterator for random numbers.
• bootstrap_iter(g,N) — bootstrap iterator.
• svd(g) — SVD modification of correlation matrix.
• dataset.bin_data(data) — bin random sample data.
• dataset.avg_data(data) — estimate means of random sample data.
• dataset.bootstrap_iter(data,N) — bootstrap random sample data.
• class dataset.Dataset — class for collecting random sample data.

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 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.

Args::

g (dictionary, array, or gvar.GVar): Collection of gvar.GVars whose

correlation matrix is to be determined.

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-15, fmt=False)

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 with 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)

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

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

Unlike pickle, 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.
• use_json (bool) – Data are serialized using pickle if False or json if True.

gvar.dumps(g)

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

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

Unlike pickle, 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.
• use_json (bool) – Data are serialized using pickle if False or json if True.

gvar.load(inputfile)

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.
• use_json (bool) – Data assumed serialized using pickle if False or json if True. If use_json=None (default) each of pickle and json is tried.

Returns:

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

gvar.loads(inputstring)

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().
• use_json (bool) – Data assumed serialized using pickle if False or json if True. If use_json=None (default) each of pickle and json is tried.

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:

• data (BufferDict, array) – Disassembled collection of gvar.GVars from gvar.disassemble() that are to be reassembled.
• cov (gvar.smat) – Covariance matrix corresponding to the gvar.GVars in data. (Default is gvar.gvar.cov.)

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 function creates an iterator that generates random arrays 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=None)

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 (array or dictionary or BufferDict or GVar) – An array (or dictionary) of objects of type gvar.GVar; or a 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 None.

Returns:

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

gvar.bootstrap_iter(g, n=None, svdcut=None)

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 None.

Returns:

An iterator that returns bootstrap copies of g.

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-15, wgts=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.

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.

The input parameters are :

Parameters:

• g – An array of gvar.GVars or a dicitionary whose values are gvar.GVars and/or arrays of gvar.GVars.
• svdcut (None or number (|svdcut|<=1).) – 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. Default is 1e-15.
• 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.

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 of g’s covariance matrix is stored in svd itself:

svd.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.

svd.nmod

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

svd.nblocks

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

svd.eigen_range

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

svd.logdet

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

svd.correction

Array containing the svd corrections that were added to g.flat to create the modified gs.

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] (gvar.GVar or array/dictionary of gvar.GVars) – Variables 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

The following class is a specialized form of an ordered dictionary for holding gvar.GVars (or other scalars) and arrays of gvar.GVars (or other scalars) that supports Python pickling:

class gvar.BufferDict

Ordered dictionary whose data are packed into a 1-d buffer (numpy.array).

A gvar.BufferDict object is an ordered dictionary whose values must either be scalars or arrays (like numpy arrays, with arbitrary shapes). The scalars and arrays are assembled into different parts of a single one-dimensional buffer. The various scalars and arrays are retrieved using keys: e.g.,

>>> a = BufferDict()
>>> a['scalar'] = 0.0
>>> a['vector'] = [1.,2.]
>>> a['tensor'] = [[3.,4.],[5.,6.]]
>>> print(a.flatten())              # print a's buffer
[ 0.  1.  2.  3.  4.  5.  6.]
>>> for k in a:                     # iterate over keys in a
...     print(k,a[k])
scalar 0.0
vector [ 1.  2.]
tensor [[ 3.  4.]
 [ 5.  6.]]
>>> a['vector'] = a['vector']*10    # change the 'vector' part of a
>>> print(a.flatten())
[  0.  10.  20.   3.   4.   5.   6.]

The first four lines here could have been collapsed to one statement:

a = BufferDict(scalar=0.0,vector=[1.,2.],tensor=[[3.,4.],[5.,6.]])

or

a = BufferDict([('scalar',0.0),('vector',[1.,2.]),
                ('tensor',[[3.,4.],[5.,6.]])])

where in the second case the order of the keys is preserved in a (since BufferDict is an ordered dictionary).

The keys and associated shapes in a gvar.BufferDict can be transferred to a different buffer, creating a new gvar.BufferDict: e.g., using a from above,

>>> buf = numpy.array([0.,10.,20.,30.,40.,50.,60.])
>>> b = BufferDict(a, buf=buf)          # clone a but with new buffer
>>> print(b['tensor'])
[[ 30.  40.]
 [ 50.  60.]]
>>> b['scalar'] += 1
>>> print(buf)
[  1.  10.  20.  30.  40.  50.  60.]

Note how b references buf and can modify it. One can also replace the buffer in the original gvar.BufferDict using, for example, a.buf = buf:

>>> a.buf = buf
>>> print(a['tensor'])
[[ 30.  40.]
 [ 50.  60.]]
>>> a['tensor'] *= 10.
>>> print(buf)
[  1.  10.  20.  300.  400.  500.  600.]

a.buf is the numpy array used for a’s buffer. It can be used to access and change the buffer directly. In a.buf = buf, the new buffer buf must be a numpy array of the correct shape. The buffer can also be accessed through iterator a.flat (in analogy with numpy arrays), and through a.flatten() which returns a copy of the buffer.

When creating a gvar.BufferDict from a dictionary (or another gvar.BufferDict), the keys included and their order can be specified using a list of keys: for example,

>>> d = dict(a=0.0,b=[1.,2.],c=[[3.,4.],[5.,6.]],d=None)
>>> print(d)
{'a': 0.0, 'c': [[3.0, 4.0], [5.0, 6.0]], 'b': [1.0, 2.0], 'd': None}
>>> a = BufferDict(d, keys=['d', 'b', 'a'])
>>> for k in a:
...     print(k, a[k])
d None
b [1.0 2.0]
a 0.0

A gvar.BufferDict functions like a dictionary except: a) items cannot be deleted once inserted; b) all values must be either scalars or arrays of scalars, where the scalars can be any noniterable type that works with numpy arrays; and c) any new value assigned to an existing key must have the same size and shape as the original value.

Note that gvar.BufferDicts can be pickled and unpickled even when they store gvar.GVars (which themselves cannot be pickled separately).

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).

The main methods 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.

isscalar(k)

Return True if self[k] is scalar else False.

update(d)

Add contents of dictionary d to self.

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 or mat. Note that val[i]<=val[i+1].

vec

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

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.

vegas.PDFIntegrator and other PDF-related Objects

Expectation values using probability density functions defined by collections of gvar.GVars can be evaluated using the vegas module (for multi-dimensional integration) and class vegas.PDFIntegrator. Related classes are:

class gvar.PDF(g, svdcut=1e-15)

Probability density function (PDF) for g.

Given an array or dictionary g of gvar.GVars, pdf=PDF(g) is the probability density function for the (usually correlated) multi-dimensional Gaussian distribution defined by g. That is pdf(p) is the probability density for random sample p drawn from g. The logarithm of the PDF is obtained using pdf.logpdf(p).

Parameters:

• g – gvar.GVar or array of gvar.GVars, or dictionary of gvar.GVars or arrays of gvar.GVars.
• svdcut (non-negative float or None) – If not None, replace covariance matrix of g with a new matrix whose small eigenvalues are modified: eigenvalues smaller than svdcut times the maximum eigenvalue eig_max are replaced by svdcut*eig_max. This can ameliorate problems caused by roundoff errors when inverting the covariance matrix. It increases the uncertainty associated with the modified eigenvalues and so is conservative. Setting svdcut=None or svdcut=0 leaves the covariance matrix unchanged. Default is 1e-15.

class gvar.PDFStatistics(moments=None, histogram=None)

Compute statistical information about a distribution.

Given moments m[i] of a random variable, computes mean, standard deviation, skewness, and excess kurtosis.

Parameters:

• moments (array of floats) – moments[i] is the (i+1)-th moment. Optional argument unless histgram=None.
• histogram (tuple) – Tuple (bins,prob) where prob[i] is the probability in the bin between bins[i-1] and bins[i]. prob[0] is the probability below bins[0] and prob[-1] is the probability above bins[-1]. Array bins is ordered. The format for prob is what is returned by accumulating calls to gvar.PDFHistogram.count(). Optional argument unless moments=None.

The attributes are as follows:

mean

mean value

sdev

standard deviation

skew

skewness coefficient

ex_kurt

excess kurtosis

median

median (if histogram provided)

plus

interval (median, median+plus) contains 34.1% of probability

minus

interval (median-minus, median) contains 34.1% of probability

gvar

gvar.gvar(mean, sdev)

class gvar.PDFHistogram(g, nbin=None, binwidth=None, bins=None)

Utility class for creating PDF histograms.

This class is designed to facilitate studies of probability density functions associated with gvar.GVars. The following code, for example, makes a histogram of probabilities for the Gaussian distribution corresponding to gvar.GVar 1.0(5):

g = gv.gvar('1.0(5)')
data = [g() for i in range(10000)]

hist = gv.PDFHistogram(g)
count = hist.count(data)
a = hist.analyze(count)
print('probabilities:', a.prob)
print('statistics:\n', a.stats)

Here hist defines a histogram with 8 bins, centered on g.mean, and each with width equal to g.sdev. The data in array data is a random sampling from g’s distribution. The number of data elements in each bin is determined by hist.count(data) and turned into probabilities by hist.analyze(count). The probabilities (a.prob) and a statistical analysis of the probability distribution based on the histogram (a.stats) are then printed out:

probabilities: [ 0.0017  0.0213  0.1358  0.3401  0.3418  0.1351  0.0222  0.0018]
statistics:
      mean = 1.001   sdev = 0.52334   skew = 0.0069999   ex_kurt = 0.034105
   median = 1.00141666398   plus = 0.499891542549   minus = 0.501710986504

A plot of the histogram can be created and displayed using, for example:

plt = hist.make_plot(count)
plt.xlabel('g')
plt.ylabel('probabilities')
plt.show()

vegas.PDFIntegrator can be used to create histograms for more complicated, multi-dimensional distributions: the expectation value of hist.count(f(p)) over values of p drawn from a multi-dimensional distribution gives the the probability distribution for function f(p).

Parameters:

• g (gvar.GVar or None) – The mean and standard deviation of g are used to design the histogram bins, which are centered on g.mean. Ignored if None (in which case bins must be specified).
• nbin (int) – The number of histogram bins. Set equal to PDFHistogram.default_nbin (=8 initially) if None.
• binwidth (float) – The width of each bin is binwidth * g.sdev. Set equal to PDFHistogram.default_binwidth (=1. initially) if None.
• bins (array or None) – Ignored if None (default). Otherwise specifies the histogram’s bin edges, overriding the default bin design specified by g. len(bins) is one larger than the number of bins. If specified it overrides the default bin design indicated by g. One of g or bins must be specified.

The main attributes are:

g

gvar.GVar used to design the histogram.

bins

Bin edges for the histogram (see above).

midpoints

Bin midpoints.

widths

Bin widths.

The main methods are:

count(data)

Compute histogram of data.

Counts the number of elements from array data in each bin of the histogram. Results are returned in an array, call it h, of length nbin+2 where h[0] is the number of data elements that fall below the range of the histogram, h[-1] (i.e., h[nbin+1]) is the number that fall above the range, and h[i] is the number in the i-th bin for i=1...nbin.

Argument data can also be a float, in which case the result is the same as from histogram([data]). Note that the expectation value of count(f(p)) over parameter values p drawn from a random distribution gives the probabilities for values of f(p) to fall in each histogram bin. Dividing by the bin widths gives the average probability density for random variable f(p) in each bin.

Bin intervals are closed on the left and open on the right, except for the last interval which is closed on both ends.

analyze(count)

Analyze count data from PDFHistogram.count().

Turns an array of counts (see PDFHistogram.count()) into a histogram of probabilities, and estimates the mean, standard deviation, and other statistical characteristics of the corresponding probability distribution.

Parameters:count (array) – Array of length nbin+2 containing histogram data where count[0] is the count for values that are below the range of the histogram, count[-1] is the count for values above the range, and count[i] is the count for the i-th bin where i=1...nbin.

Returns a named tuple containing the following information (in order):

bins: Array of bin edges for histogram (length nbin+1)

prob: Array of probabilities for each bin.

stats: Statistical data about histogram. See PDFStatistics.

norm: Convert counts into probabilities by dividing by norm.

make_plot(count, plot=None, show=False, plottype=’probability’, bar=dict(alpha=0.15, color=’b’), errorbar=dict(fmt=’b.’), gaussian=dict(ls=’–’, c=’r’))

Convert histogram counts in array count into a plot.

Parameters:

• count (array) – Array of histogram counts (see PDFHistogram.count()).
• plot (plotter) – matplotlib plotting window. If None uses the default window. Default is None.
• show (boolean) – Displayes plot if True; otherwise returns the plot. Default is False.
• plottype (str) – The probabilities in each bin are plotted if plottype='probability' (default). The average probability density is plot if plottype='density'. The cumulative probability is plotted if plottype=cumulative.
• bar (dictionary) – Additional plotting arguments for the bar graph showing the histogram. This part of the plot is omitted if bar=None.
• errorbar (dictionary) – Additional plotting arguments for the errorbar graph, showing error bars on the histogram. This part of the plot is omitted if errorbar=None.
• gaussian (dictionary) – Additional plotting arguments for the plot of the Gaussian probability for the gvar.GVar (g) specified in the initialization. This part of the plot is omitted if gaussian=None or if no g was specified.

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.