Python documentation

The package chiexp can be loaded in python using the standard import command:

1import sys
2sys.path.append('/path/to/chiexp/directory/lib/python')
3from chiexp import chisquare

The chisquare class

class chiexp.chisquare(x, y, W, f, df, v='x')

A class with the relevant functionalities to compute the expected chi square.

Parameters:

  • x (array) – array with the values of the x-axis; 2-D arrays are accepted and the second dimension is intepreted as internal kinematic index

  • y (array) – array with the values of the y-axis, of same length as x

  • W (array) – the weight matrix, N-by-N, or its diagional of length N

  • f (function) – callable function or lambda function defining the fitted function; the program assumes x correspond to the first arguments

  • df (function) – callable function or lambda function returning an array that contains the gradient of f, namely \(\partial \phi(\{p\},\{x_i\})/\partial p_\alpha\)

  • v (str, optional) – a string with the list of variables used in f as the kinematic coordinates. Default value corresponds to x, which implies that f must be defined using x as first and unique kinematic variable.

chiexp(cov, Nrep=None, Stau=1.5, Wopt=None, Wcov=None, plot=False)

Computes the expected chi square

Parameters:

  • cov (list or array) – the covariance matrix is assumed if the object passed has is a N-by-N 2D array; otherwise if it is a M-by-N 2D array the program assumes that it contains the fluctuations of the N observables over M configurations

  • Nrep (list, optional) – number of configurations per replica, e.g. if Nrep=[N1,N2,N3] then M=N1+N2+N3

  • Stau (float, optional) – parameter used in the automatic window procedure

  • Wopt (float, optional) – optimal window used in the estimate of expected chi square; if passed Stau is ignored

  • Wcov (float, optional) – window used in the estimate of all entries of the covariance matrix

  • plot (bool, optional) – if set to True the program produces a plot with the autocorrelation function of the expected chi square

Returns:

the expected chi square, its error and the estimated covariance matrix; if Wcov is not

passed the window obtained from the autocorrelation function of the expected chi square is used for all elements; if cov is a N-by-N 2D array it returns a copy of it

Return type:

list

Note

The additional parameters Nrep,`Stau` and plot are ignored if cov corresponds to the covariance matrix, i.e. it is a N-by-N array. The matrix \([C^{1/2} W^{1/2} (1-P) W^{1/2} C^{1/2}]\) can be accessed by invoking the nu element of the chisquare class

chisq(p)

Returns the chi square from the input parameters p.

derfit()

Returns the derivative of the parameters w.r.t. Y

At the minimum the parameters depend on the input observables; their errors can be obtained by error propagation through the derivatives provided by this routine.

Returns:

the derivatives der(i,a) = dp(a)/dy(i)

Return type:

array

fit(p0, min_search)

Minimizes the chi square starting from the initial guess parameters p0

Parameters:

  • p0 (array) – the initial guess values of the parameters

  • min_search (function) – an external minimizer

Returns:

output parameters and the value of the chi square at the minimum

Return type:

list

pvalue(method='eig', nmc=5000, plot=False)

Computes the p-value of the fit

Parameters:

  • method (string, optional) – string specifying the method to estimate the quality of fit. Accepted values are ‘MC’ for a pure Monte Carlo estimate or ‘eig’ (default) for the formula based on the eigenvalues of the matrix nu.

  • nmc (int, optional) – number of Monte Carlo samples used to estimate the quality of fit. Default is 5000.

  • plot (bool, optional) – if set to True plots the probabilty distribution of the expected chi square

Returns:

the quality of fit, the error of the quality of fit and the Monte Carlo history of the expected chi square

Return type:

list

Note

The class must know the value of the chi square at the minimum, which means that either fit or chisq must be called before qfit. If method is set to ‘MC’ the error of the quality of fit is based only the MC sampling.

set_pars(p0)

Sets the parameters to the input values p0

Examples of usage

First we define a new instance of the chisquare class; in the following example we do it for an uncorrelated fit

1[x, y, dy] = load_your_data_set()
2
3# if dy is the error of y we define W for an uncorrelated fit
4W = [1./e**2 for e in dy]
5
6func = lambda x, a, m: a*exp(-m*x)
7dfunc = lambda x, a, m: [exp(-m*x), -a*x*exp(-m*x)]
8
9c=chisquare(x,y,W,func,dfunc)

If the minimum of \(\chi^2\) is not know the user can use the fit method; otherwise the values of the parameters at the minimum can be passed to the class via the chisq method

1>>> from scipy.optimize import minimize
2>>> p0=[1.,1.] # guess values of the paramters
3>>> [p, c2] = c.fit(p0, minimize)
4>>> # or pass the parameters at the minimum
5>>> c.chisq(p)

The user can also compute the error of the fitted parameters from the errors (fluctuations) of the input observables \(Y\). To do so the derivatives \(d p_\alpha / d y_i\) (defined at the minimum of the \(\chi^2\)) are needed and if the covariance matrix is known, applying the chain rule returns the wanted errors

>>> der = c.derfit() # N x Na matrix
>>> print (der.T @ cov @ der) # errors of parameters

\(\langle \chi^2 \rangle\) can be computed using the method chiexp

1>>> [ce, dce, _] = c.chiexp(cov)
2>>> # if cov contains the fluctuations of M configs and 2 replica
3>>> (M, N) = numpy.shape(cov)
4>>> nr = [N1 N2] # such that M=N1+N2
5>>> print M - sum(nr)
60
7>>> [ce,dce,covest] = c.chiexp(cov,Nrep=nr,Stau=2.5)

Finally the quality of fit is estimated from a Monte-Carlo chain of length nmc

>>> nmc=10000
>>> [p, dp, h] = c.pvalue(nmc, plot=True)

The method pvalue returns the quality of fit and its error, p and dp respectively, and the Monte Carlo history of \(\chi^2\) h. If the plot flag is activated a plot is automatically generated with the distribution probability, namely a normalized histogram of h.