lsqfit
- Nonlinear Least Squares Fitting¶
Introduction¶
This package contains tools for nonlinear least-squares curve fitting of data. In general a fit has four inputs:
The dependent data
y
that is to be fit — typicallyy
is a Python dictionary in anlsqfit
analysis. Its valuesy[k]
are eithergvar.GVar
s or arrays (any shape or dimension) ofgvar.GVar
s that specify the values of the dependent variables and their errors.A collection
x
of independent data —x
can have any structure and contain any data, or it can be omitted.A fit function
f(x, p)
whose parametersp
are adjusted by the fit untilf(x, p)
equalsy
to withiny
s errors — parameters p` are usually specified by a dictionary whose valuesp[k]
are individual parameters or (numpy
) arrays of parameters. The fit function is assumed independent ofx
(that is,f(p)
) ifx = False
(or ifx
is omitted from the input data).Initial estimates or priors for each parameter in
p
— priors are usually specified using a dictionaryprior
whose valuesprior[k]
aregvar.GVar
s or arrays ofgvar.GVar
s that give initial estimates (values and errors) for parametersp[k]
.
A typical code sequence has the structure:
... collect x, y, prior ...
def f(x, p):
... compute fit to y[k], for all k in y, using x, p ...
... return dictionary containing the fit values for the y[k]s ...
fit = lsqfit.nonlinear_fit(data=(x, y), prior=prior, fcn=f)
print(fit) # variable fit is of type nonlinear_fit
The parameters p[k]
are varied until the chi**2
for the fit is
minimized.
The best-fit values for the parameters are recovered after fitting
using, for example, p=fit.p
. Then the p[k]
are gvar.GVar
s or
arrays of gvar.GVar
s that give best-fit estimates and fit uncertainties
in those estimates (as well as the correlations between them).
The print(fit)
statement prints a summary of the fit results.
The dependent variable y
above could be an array instead of a
dictionary, which is less flexible in general but possibly more
convenient in simpler fits. Then the approximate y
returned by fit
function f(x, p)
must be an array with the same shape as the dependent
variable. The prior prior
could also be represented by an array
instead of a dictionary.
By default priors are Gaussian/normal distributions, represented by
gvar.GVar
s. lsqfit
also
allows for log-normal and other distributions as well. The
latter are indicated by replacing the prior (in a dictionary prior)
with key c
, for example, by a prior for the parameter’s logarithm,
with key log(c)
.
nonlinear_fit
in effect adds parameter c
to the parameter
dictionary, deriving its value from parameter log(c)
.
The fit function can be expressed directly in terms of
parameter c
and so is the same no matter which distribution is
used for c
. Additional distributions
can be added using gvar.BufferDict.add_distribution()
.
The lsqfit
tutorial contains extended explanations and examples.
The first appendix in the paper at http://arxiv.org/abs/arXiv:1406.2279
provides conceptual background on the techniques used in this
module for fits and, especially, error budgets.
nonlinear_fit Objects¶
- class lsqfit.nonlinear_fit(data, fcn, prior=None, p0=None, svdcut=1e-12, eps=None, noise=False, debug=False, tol=1e-8, maxit=1000, fitter='gsl_multifit', **fitterargs)¶
Nonlinear least-squares fit.
lsqfit.nonlinear_fit
fits a (nonlinear) functionf(x, p)
to datay
by varying parametersp
, and stores the results: for example,fit = nonlinear_fit(data=(x, y), fcn=f, prior=prior) # do fit print(fit) # print fit results
The best-fit values for the parameters are in
fit.p
, while thechi**2
, the number of degrees of freedom, the logarithm of Gaussian Bayes Factor, the number of iterations (or function evaluations), and the cpu time needed for the fit are infit.chi2
,fit.dof
,fit.logGBF
,fit.nit
, andfit.time
, respectively. Results for individual parameters infit.p
are of typegvar.GVar
, and therefore carry information about errors and correlations with other parameters. The fit data and prior can be recovered usingfit.x
(equalsFalse
if there is nox
),fit.y
, andfit.prior
; the data and prior are corrected for the SVD cut, if there is one (that is, their covariance matrices have been modified in accordance with the SVD cut).- Parameters:
data (dict, array or tuple) –
Data to be fit by
lsqfit.nonlinear_fit
can have one of the following forms:data = x, y
x
is the independent data that is passed to the fit function with the fit parameters:fcn(x, p)
.y
is a dictionary (or array) ofgvar.GVar
s that encode the means and covariance matrix for the data that is to be fit being fit. The fit function must return a result having the same layout asy
.data = y
y
is a dictionary (or array) ofgvar.GVar
s that encode the means and covariance matrix for the data being fit. There is no independent data so the fit function depends only upon the fit parameters:fit(p)
. The fit function must return a result having the same layout asy
.
Setting
x=False
in the first of these formats implies that the fit function depends only on the fit parameters: that is,fcn(p)
instead offcn(x, p)
. (This is not assumed ifx=None
.)fcn (callable) – The function to be fit to
data
. It is either a function of the independent datax
and the fit parametersp
(fcn(x, p)
), or a function of just the fit parameters (fcn(p)
) when there is nox
data orx=False
. The parameters are tuned in the fit until the function returns values that agree with they
data to within they
s’ errors. The function’s return value must have the same layout as they
data (a dictionary or an array). The fit parametersp
are either: 1) a dictionary where eachp[k]
is a single parameter or an array of parameters (any shape); or, 2) a single array of parameters. The layout of the parameters is the same as that of priorprior
if it is specified; otherwise, it is inferred from of the starting valuep0
for the fit.prior (dict, array, str, gvar.GVar or None) – A dictionary (or array) containing a priori estimates for all parameters
p
used by fit functionfcn(x, p)
(orfcn(p)
). Fit parametersp
are stored in a dictionary (or array) with the same keys and structure (or shape) asprior
. The default value isNone
;prior
must be defined ifp0
isNone
.p0 (dict, array, float, None, or True) – Starting values for fit parameters in fit.
lsqfit.nonlinear_fit
adjustsp0
to make it consistent in shape and structure withprior
when the latter is specified: elements missing fromp0
are filled in usingprior
, and elements inp0
that are not inprior
are discarded. Ifp0
is a string, it is taken as a file name andlsqfit.nonlinear_fit
attempts to read starting values from that file; best-fit parameter values are written out to the same file after the fit (for priming future fits). Ifp0
isNone
or the attempt to read the file fails, starting values are extracted fromprior
. Ifp0
isTrue
, it is replaced by a starting point drawn at random from theprior
distribution. The default value isNone
;p0
must be explicitly specified ifprior
isNone
.linear (list or None) – Optional list of fit parameters that appear linearly in the fit function. The fit function can be reexpressed (using variable projection) as a function that is independent of its linear parameters. The resulting fit has fewer fit parameters and typically will converge in fewer iterations, but each iteration will take longer. Whether or not the fit is faster or more robust in any particular application is a matter for experiment, but answers should be the same either way. The linear parameters are reconstructed from the nonlinear parameters (and the data) after the fit. Parameter
linear
is either: a list of dictionary keys corresponding to linear parameters when the parameters are stored in a dictionary (seeprior
); or, a list of indices corresponding to these parameters when they are stored in an array. Note that this feature is experimental; the interface may change in the future.svdcut (float) – If nonzero, singularities in the correlation matrix for
y
andprior
are regulated usinggvar.regulate()
with an SVD cutoffsvdcut
. This makes the correlation matrices less singular, which can improve the stability and accuracy of a fit. Default issvdcut=1e-12
.eps (float) – If positive, singularities in the correlation matrix for
y
andprior
are regulated usinggvar.regulate()
with cutoffeps
. This makes the correlation matrices less singular, which can improve the stability and accuracy of a fit. Ignored ifsvdcut
is specified (and notNone
).noise (tuple or bool) – If
noise[0]=True
, noise is added to the data means commensurate with the additional uncertainties introduced by usingsvdcut>0
oreps>0
. Ifnoise[1]=True
, noise is added to the prior means commensurate with the uncertainties in the prior. Noise is useful for testing the quality of a fit (chi2
). Settingnoise=True
is shorthand fornoise=(True, True)
, andnoise=False
meansnoise=(False, False)
(the default).udata (dict, array or tuple) – Same as
data
but instructs the fitter to ignore correlations between different pieces of data. This speeds up the fit, particularly for large amounts of data, but ignores potentially valuable information if the data actually are correlated. Only one ofdata
orudata
should be specified. (Default isNone
.)fitter (str or None) – Fitter code. Options if GSL is installed include:
'gsl_multifit'
(default) and'gsl_v1_multifit'
(original fitter). Options ifscipy
is installed include:'scipy_least_squares'
(default if GSL not installed).gsl_multifit
has many options, providing extensive user control.scipy_least_squares
can be used for fits where the parameters are bounded. (Bounded parameters can also be implemented, for any of the fitters, using non-Gaussian priors — see the tutorial.)tol (float or tuple) –
Assigning
tol=(xtol, gtol, ftol)
causes the fit to stop searching for a minimum when any ofxtol >=
relative change in parameters between iterationsgtol >=
relative size of gradient ofchi**2
functionftol >=
relative change inchi**2
between iterations
is satisfied. See the fitter documentation for detailed definitions of these stopping conditions. Typically one sets
xtol=1/10**d
whered
is the number of digits of precision desired in the result, whilegtol<<1
andftol<<1
. Settingtol=delta
wheredelta
is a number is equivalent to settingtol=(delta,1e-10,1e-10)
. Settingtol=(delta1,delta2)
is equivalent to settingtol=(delta1,delta2,1e-10)
. Default istol=1e-8
. (Note: theftol
option is disabled in some versions of the GSL library.)maxit (int) – Maximum number of algorithm iterations (or function evaluations for some fitters) in search for minimum; default is 1000.
debug (bool) – Set to
True
for extra debugging of the fit function and a check for roundoff errors. (Default isFalse
.)fitterargs (dict) – Dictionary of additional arguments passed through to the underlying fitter. Different fitters offer different parameters; see the documentation for each.
Objects of type
lsqfit.nonlinear_fit
have the following attributes:- chi2¶
The minimum
chi**2
for the fit.fit.chi2 / fit.dof
is usually of order one in good fits. Values much less than one suggest that actual fluctuations in the input data and/or priors might be smaller than suggested by the standard deviations (or covariances) used in the fit.- Type:
float
- cov¶
Covariance matrix of the best-fit parameters from the fit.
- Type:
array
- dof¶
Number of degrees of freedom in the fit, which equals the number of pieces of data being fit when priors are specified for the fit parameters. Without priors, it is the number of pieces of data minus the number of fit parameters.
- Type:
int
- error¶
Error message generated by the underlying fitter when an error occurs.
None
otherwise.- Type:
str
- fitter_results¶
Results returned by the underlying fitter. Refer to the appropriate fitter’s documentation for details.
- logGBF¶
The logarithm of the probability (density) of obtaining the fit data by randomly sampling the parameter model (priors plus fit function) used in the fit — that is, it is
P(data|model)
. This quantity is useful for comparing fits of the same data to different models, with different priors and/or fit functions. The model with the largest value offit.logGBF
is the one preferred by the data. The exponential of the difference infit.logGBF
between two models is the ratio of probabilities (Bayes factor) for those models. Differences infit.logGBF
smaller than 1 are not very significant. Gaussian statistics are assumed when computingfit.logGBF
.- Type:
float or None
- p¶
Best-fit parameters from fit. Depending upon what was used for the prior (or
p0
), it is either: a dictionary (gvar.BufferDict
) ofgvar.GVar
s and/or arrays ofgvar.GVar
s; or an array (numpy.ndarray
) ofgvar.GVar
s.fit.p
represents a multi-dimensional Gaussian distribution which, in Bayesian terminology, is the posterior probability distribution of the fit parameters.- Type:
dict, array or gvar.GVar
- pmean¶
Means of the best-fit parameters from fit.
- Type:
dict, array or float
- psdev¶
Standard deviations of the best-fit parameters from fit.
- Type:
dict, array or float
- palt¶
Same as
fit.p
except that the errors are computed directly fromfit.cov
. This is faster but means that no information about correlations with the input data is retained (unlike infit.p
); and, therefore,fit.palt
cannot be used to generate error budgets.fit.p
andfit.palt
give the same means and normally give the same errors for each parameter. They differ only when the input data’s covariance matrix is too singular to invert accurately (because of roundoff error), in which case refitting with a nonzero value foreps
orsvdcut
is advisable.- Type:
dict, array or gvar.GVar
- p0¶
The parameter values used to start the fit. This will differ from the input
p0
if the latter was incomplete.- Type:
dict, array or float
- prior¶
Prior used in the fit. This may differ from the input prior if an SVD cut is used. It is either a dictionary (
gvar.BufferDict
) or an array (numpy.ndarray
), depending upon the input. EqualsNone
if no prior was specified.- Type:
dict, array, gvar.GVar or None
- Q¶
The probability that the
chi**2
from the fit could have been larger, by chance, assuming the best-fit model is correct. Good fits haveQ
values larger than 0.1 or so. Also called the p-value of the fit. The probabilistic intrepretation becomes unreliable if the actual fluctuations in the input data and/or priors are much smaller than suggested by the standard deviations (or covariances) used in the fit (leading to an unusually smallchi**2
).- Type:
float or None
- residuals¶
An array containing the fit residuals normalized by the corresponding standard deviations. The residuals are projected onto the eigenvectors of the correlation matrix and so should be uncorrelated from each other. The residuals include contributions from both the fit data and the prior. They are related to the the
chi**2
of the fit by:chi2 = sum(fit.residuals**2)
.
- stopping_criterion¶
Criterion used to stop fit:
0: didn’t converge
1:
xtol >=
relative change in parameters between iterations2:
gtol >=
relative size of gradient ofchi**2
3:
ftol >=
relative change inchi**2
between iterations4: unable to improve fit further (e.g., already converged)
- Type:
int
- correction¶
Sum of all corrections, if any, added to the fit data and prior when
eps>0
orsvdcut>0
.- Type:
gvar.GVar
- svdn¶
Number of eigenmodes of the correlation matrix modified (and/or deleted) when
svdcut>0
.- Type:
int
- time¶
CPU time (in secs) taken by fit.
- Type:
float
- tol¶
Tolerance used in fit. This differs from the input tolerance if the latter was incompletely specified.
- Type:
tuple
- x¶
The first field in the input
data
. This is sometimes the independent variable (as in ‘y vs x’ plot), but may be anything. It is set equal toFalse
if thex
field is omitted from the inputdata
. (This also means that the fit function has nox
argument: sof(p)
rather thanf(x,p)
.)- Type:
obj
- y¶
Fit data used in the fit. This may differ from the input data if an SVD cut is used. It is either a dictionary (
gvar.BufferDict
) or an array (numpy.ndarray
), depending upon the input.- Type:
dict, array or gvar.GVar
- nblocks¶
nblocks[s]
equals the number of block-diagonal sub-matrices of they
–prior
covariance matrix that are sizes
-by-s
. This is sometimes useful for debugging.- Type:
dict
The global defaults used by
lsqfit.nonlinear_fit
can be changed by changing entries in dictionarylsqfit.nonlinear_fit.DEFAULTS
for keys'eps'
,'svdcut'
,'debug'
,'tol'
,'noise'
,'maxit'
, and'fitter'
. Additional defaults can be added to that dictionary to be are passed throughlsqfit.nonlinear_fit
to the underlying fitter (via dictionaryfitterargs
).Additional methods are provided for printing out detailed information about the fit, evaluating
chi**2
, testing fits with simulated data, doing bootstrap analyses of the fit errors, dumping (for later use) and loading parameter values, and checking for roundoff errors in the final error estimates:- format(maxline=0, pstyle='v')¶
Formats fit output details into a string for printing.
The output tabulates the
chi**2
per degree of freedom of the fit (chi2/dof
), the number of degrees of freedom, theQ
value of the fit (ie, p-value), and the logarithm of the Gaussian Bayes Factor for the fit (logGBF
). At the end it lists the SVD cut, the number of eigenmodes modified by the SVD cut, the tolerances used in the fit, and the time in seconds needed to do the fit. The tolerance used to terminate the fit is marked with an asterisk. It also lists information about the fitter used if it is other than the standard choice.Optionally,
format
will also list the best-fit values for the fit parameters together with the prior for each (in[]
on each line). Lines for parameters that deviate from their prior by more than one (prior) standard deviation are marked with asterisks, with the number of asterisks equal to the number of standard deviations (up to five). Lines for parameters designated as linear (seelinear
keyword) are marked with a minus sign after their prior.format
can also list all of the data and the corresponding values from the fit, again with asterisks on lines where there is a significant discrepancy.- Parameters:
maxline (int or bool) – Maximum number of data points for which fit results and input data are tabulated.
maxline<0
implies that onlychi2
,Q
,logGBF
, anditns
are tabulated; no parameter values are included. Settingmaxline=True
prints all data points; setting it equalNone
orFalse
is the same as setting it equal to-1
. Default ismaxline=0
.pstyle (str or None) – Style used for parameter list. Supported values are ‘vv’ for very verbose, ‘v’ for verbose, and ‘m’ for minimal. When ‘m’ is set, only parameters whose values differ from their prior values are listed. Setting
pstyle=None
implies no parameters are listed.extend (bool) – If
True
, extend the parameter list to include values derived from log-normal or other non-Gaussian parameters. So values for fit parameterp['log(a)']
, for example, are listed together with valuesp['a']
for the exponential of the fit parameter. Settingextend=False
means that only the value forp['log(a)']
is listed. Default isTrue
.
- Returns:
String containing detailed information about fit.
- dchi2(p)¶
chi**2(p) - fit.chi2
for fit parametersp
.- Paramters:
- p: Array or dictionary containing values for fit parameters, using
the same layout as in the fit function.
- Returns:
chi**2(p) - fit.chi2
wherechi**2(p)
is the fit’schi**2
for fit parametersp
andfit.chi2
is thechi**2
value for the best fit.
- pdf(p)¶
exp(-(chi**2(p) - fit.chi2)/2)
for fit parametersp
.fit.pdf(p)
is proportional to the probability density function (PDF) used in the fit:fit.pdf(p)/exp(fit.pdf.lognorm)
is the product of the Gaussian PDF for the dataP(data|p,M)
times the Gaussian PDF for the priorP(p|M)
whereM
is the model used in the fit (i.e., the fit function and prior). The product of PDFs isP(data,p|M)
by Bayes’ Theorem; integrating over fit parameters p gives the Bayes Factor or EvidenceP(data|M)
, which is proportional to the probability that the fit data come from fit modelM
. The logarithm of the Bayes Factor should agree withfit.logGBF
when the Gaussian approximation assumed in the fit is accurate.fit.pdf(p)
is useful for checking a least-squares fit against the corresponding Bayesian integrals. In the following example,vegas.PDFIntegrator
from thevegas
module is used to evaluate Bayesian expectation values ofs*g
and its standard deviation wheres
andg
are fit parameters:import gvar as gv import lsqfit import numpy as np import vegas def main(): # least-squares fit x = np.array([0.1, 1.2, 1.9, 3.5]) y = gv.gvar(['1.2(1.0)', '2.4(1)', '2.0(1.2)', '5.2(3.2)']) prior = gv.gvar(dict(a='0(5)', s='0(2)', g='2(2)')) fit = lsqfit.nonlinear_fit(data=(x,y), prior=prior, fcn=fitfcn, debug=True) print(fit) # create integrator and adapt it to PDF (warmup) neval = 10_000 nitn = 10 expval = vegas.PDFIntegrator(fit.p, pdf=fit.pdf, nproc=4) warmup = expval(neval=neval, nitn=nitn) # calculate expectation value of g(p) results = expval(g, neval=neval, nitn=nitn, adapt=False) print(results.summary(True)) print('results =', results, '\n') sg = results['sg'] sg2 = results['sg2'] sg_sdev = (sg2 - sg**2) ** 0.5 print('s*g from Bayes integral: mean =', sg, ' sdev =', sg_sdev) print('s*g from fit:', fit.p['s'] * fit.p['g']) print() print('logBF =', np.log(results.pdfnorm) - fit.pdf.lognorm) def fitfcn(x, p): return p['a'] + p['s'] * x ** p['g'] def g(p): sg = p['s'] * p['g'] return dict(sg=sg, sg2=sg**2) if __name__ == '__main__': main()
Here the probability density function used for the expectation values is
fit.pdf(p)
, and the expectation values are returned in dictionaryresults
.vegas
uses adaptive Monte Carlo integration. Thewarmup
calls to the integrator are used to adapt it to the probability density function, and then the adapted integrator is called again to evaluate the expectation value. Parameterneval
is the (approximate) number of function calls per iteration of thevegas
algorithm andnitn
is the number of iterations. We use the integrator to calculated the expectation value ofs*g
and(s*g)**2
so we can compute a mean and standard deviation.The output from this code shows that the Gaussian approximation for
s*g
(0.78(66)) is somewhat different from the result obtained from a Bayesian integral (0.49(53)):Least Square Fit: chi2/dof [dof] = 0.32 [4] Q = 0.87 logGBF = -9.2027 Parameters: a 1.61 (90) [ 0.0 (5.0) ] s 0.62 (81) [ 0.0 (2.0) ] g 1.2 (1.1) [ 2.0 (2.0) ] Settings: svdcut/n = 1e-12/0 tol = (1e-08*,1e-10,1e-10) (itns/time = 18/0.0) itn integral average chi2/dof Q ------------------------------------------------------- 1 0.954(11) 0.954(11) 0.00 1.00 2 0.9708(99) 0.9622(74) 0.74 0.53 3 0.964(12) 0.9627(63) 0.93 0.47 4 0.9620(93) 0.9626(52) 0.86 0.56 5 0.964(14) 0.9629(50) 0.71 0.74 6 0.957(17) 0.9619(50) 0.65 0.84 7 0.964(12) 0.9622(46) 0.61 0.90 8 0.9367(86) 0.9590(42) 0.80 0.73 9 0.9592(94) 0.9591(39) 0.75 0.80 10 0.952(13) 0.9584(37) 0.72 0.85 key/index value ------------------------------------ pdf 0.9584 (37) ('f(p)*pdf', 'sg') 0.4652 (23) ('f(p)*pdf', 'sg2') 0.5073 (33) results = {'sg': 0.4854(20), 'sg2': 0.5293(33)} s*g from Bayes integral: mean = 0.4854(20) sdev = 0.5420(17) s*g from fit: 0.78(66) logBF = -9.1505(39)
The result
logBF
for the logarithm of the Bayes Factor from the integral agrees well withfit.logGBF
, the log Bayes Factor in the Gaussian approximation. This is evidence that the Gaussian approximation implicit in the least squares fit is reliable; the product ofs*g
, however, is not so Gaussian because of the large uncertainties (compared to the means) ins
andg
separately.- Paramters:
- p: Array or dictionary containing values for fit parameters, using
the same layout as in the fit function.
- Returns:
exp(-(chi**2(p) - fit.chi2)/2)
wherechi**2(p)
is the fit’schi**2
for fit parametersp
andfit.chi2
is thechi**2
value for the best fit.
- simulated_fit_iter(n=None, pexact=None, add_priornoise=False, **kargs)¶
Iterator that returns simulation copies of a fit.
Fit reliability is tested using simulated data which replaces the mean values in
self.y
with random numbers drawn from a distribution whose mean equalsself.fcn(pexact)
and whose covariance matrix is the same asself.y
’s. Simulated data is very similar to the original fit data,self.y
, but corresponds to a world where the correct values for the parameters (i.e., averaged over many simulated data sets) are given bypexact
.pexact
is usually taken equal tofit.pmean
.Each iteration of the iterator creates new simulated data, with different random numbers, and fits it, returning the the
lsqfit.nonlinear_fit
that results. The simulated data has the same covariance matrix asfit.y
. Typical usage is:... fit = nonlinear_fit(...) ... for sfit in fit.simulated_fit_iter(n=3): ... verify that sfit has a good chi**2 ... ... verify that sfit.p agrees with pexact=fit.pmean within errors ...
Only a few iterations are needed to get a sense of the fit’s reliability since we know the correct answer in each case. The simulated fit’s output results should agree with
pexact
(=fit.pmean
here) within the simulated fit’s errors.Setting parameter
add_priornoise=True
varies the means of the priors as well as the means of the data. This option is useful for testing goodness of fit because with itchi**2/N
should be1 ± sqrt(2/N)
, whereN
is the number of degrees of freedom. (chi**2/N
can be significantly smaller than one without added noise in prior means.)Simulated fits can also be used to estimate biases in the fit’s output parameters or functions of them, should non-Gaussian behavior arise. This is possible, again, because we know the correct value for every parameter before we do the fit. Again only a few iterations may be needed for reliable estimates.
- Parameters:
n (int or
None
) – Maximum number of iterations (equals infinity ifNone
).pexact (
None
or array/dict of numbers) – Fit-parameter values for the underlying distribution used to generate simulated data; replaced byself.pmean
if isNone
(default).add_priornoise (bool) – Vary prior means if
True
; otherwise vary only the means inself.y
(default).kargs – Dictionary containing override values for fit parameters.
- Returns:
An iterator that returns
lsqfit.nonlinear_fit
s for different simulated data.
- simulated_data_iter(n=None, pexact=None, add_priornoise=False)¶
Iterator that returns simulated data based upon a fit’s data.
Simulated data is generated from a fit’s data
fit.y
by replacing the mean values in that data with random numbers drawn from a distribution whose mean isself.fcn(pexact)
and whose covariance matrix is the same as that ofself.y
. Each iteration of the iterator returns new simulated data, with different random numbers for the means and a covariance matrix equal to that ofself.y
. This iterator is used byself.simulated_fit_iter
.Typical usage:
fit = nonlinear_fit(data=(x,y), prior=prior, fcn=fcn) ... for ysim, priorsim in fit.simulate_data_iter(n=10): fitsim = nonlinear_fit(data=(x, ysim), prior=priorsim, fcn=fcn) print(fitsim) print('chi2 =', gv.chi2(fit.p, fitsim.p))
This code tests the fitting protocol on simulated data, comparing the best fit parameters in each case with the correct values (
fit.p
). The loop in this code is functionally the same as (but probably not as fast as):for fitsim in fit.simulated_fit_iter(n=10): print(fitsim) print('chi2 =', gv.chi2(fit.p, fitsim.p))
- Parameters:
n (int or None) – Maximum number of iterations (equals infinity if
None
).pexact (None or dict/array of numbers) – Fit-parameter values for the underlying distribution used to generate simulated data; replaced by
self.pmean
if isNone
(default).add_priornoise (bool) – Vary prior means if
True
; otherwise vary only the means inself.y
(default).
- Returns:
An iterator that returns a 2-tuple containing simulated versions of self.y and self.prior:
(ysim, priorsim)
.
- bootstrapped_fit_iter(n=None, datalist=None)¶
Iterator that returns bootstrap copies of a fit.
A bootstrap analysis involves three steps: 1) make a large number of “bootstrap copies” of the original input data and prior that differ from each other by random amounts characteristic of the underlying randomness in the original data; 2) repeat the entire fit analysis for each bootstrap copy of the data, extracting fit results from each; and 3) use the variation of the fit results from bootstrap copy to bootstrap copy to determine an approximate probability distribution (possibly non-gaussian) for the fit parameters and/or functions of them: the results from each bootstrap fit are samples from that distribution.
Bootstrap copies of the data for step 2 are provided in
datalist
. Ifdatalist
isNone
, they are generated instead from the means and covariance matrix of the fit data (assuming gaussian statistics). The maximum number of bootstrap copies considered is specified byn
(None
implies no limit).Variations in the best-fit parameters (or functions of them) from bootstrap fit to bootstrap fit define the probability distributions for those quantities. For example, one could use the following code to analyze the distribution of function
g(p)
of the fit parameters:fit = nonlinear_fit(...) ... glist = [] for bsfit in fit.bootstrapped_fit_iter( n=100, datalist=datalist, ): glist.append(g(bsfit.pmean)) ... analyze samples glist[i] from g(p) distribution ...
This code generates
n=100
samplesglist[i]
from the probability distribution ofg(p)
. If everything is Gaussian, the mean and standard deviation ofglist[i]
should agree withg(fit.p).mean
andg(fit.p).sdev
.- Parameters:
n (int) – Maximum number of iterations if
n
is notNone
; otherwise there is no maximum.datalist (iter) – Collection of bootstrap
data
sets for fitter.kargs (dict) – Overrides arguments in original fit.
- Returns:
Iterator that returns an
lsqfit.nonlinear_fit
object containing results from the fit to the next data set indatalist
.
- check_roundoff(rtol=0.25, atol=1e-6)¶
Check for roundoff errors in fit.p.
Compares standard deviations from fit.p and fit.palt to see if they agree to within relative tolerance
rtol
and absolute toleranceatol
. Generates a warning if they do not (in which case an SVD cut might be advisable).
- qqplot_residuals(plot=None)¶
QQ plot normalized fit residuals.
The sum of the squares of the residuals equals
self.chi2
. Individual residuals should be distributed in a Gaussian distribution centered about zero. A Q-Q plot orders the residuals and plots them against the value they would have if they were distributed according to a Gaussian distribution. The resulting plot will approximate a straight line along the diagonal of the plot (dashed black line) if the residuals have a Gaussian distribution with zero mean and unit standard deviation.The residuals are fit to a straight line and the fit is displayed in the plot (solid red line). Residuals that fall on a straight line have a distribution that is Gaussian. A nonzero intercept indicates a bias in the mean, away from zero. A slope smaller than 1.0 indicates the actual standard deviation is smaller than suggested by the fit errors, as would be expected if the
chi2/dof
is significantly below 1.0 (sincechi2
equals the sum of the squared residuals).One way to display the plot is with:
fit.qqplot_residuals().show()
- Parameters:
plot – a
matplotlib
plotter. IfNone
, usesmatplotlib.pyplot
.- Returns:
Plotter
plot
.
This method requires the
scipy
andmatplotlib
modules.
- plot_residuals(plot=None)¶
Plot normalized fit residuals.
The sum of the squares of the residuals equals
self.chi2
. Individual residuals should be distributed about one, in a Gaussian distribution.- Parameters:
plot –
matplotlib
plotter. IfNone
, usesmatplotlib.pyplot
.- Returns:
Plotter
plot
.
- static set(clear=False, **defaults)¶
Set default parameters for
lsqfit.nonlinear_fit
.Use to set default values for parameters:
eps
,svdcut
,debug
,tol
,maxit
, andfitter
. Can also set parameters specific to the fitter specified by thefitter
argument.Sample usage:
import lsqfit old_defaults = lsqfit.nonlinear_fit.set( fitter='gsl_multifit', alg='subspace2D', solver='cholesky', tol=1e-10, debug=True, )
nonlinear_fit.set()
without arguments returns a dictionary containing the current defaults.- Parameters:
clear (bool) – If
True
remove earlier settings, restoring the original defaults, before adding new defaults. The default value isclear=False
.nonlinear_fit.set(clear=True)
restores the original defaults.defaults (dict) – Dictionary containing new defaults.
- Returns:
A dictionary containing the old defaults, before they were updated. These can be restored using
nonlinear_fit.set(old_defaults)
whereold_defaults
is the dictionary containint the old defaults.
vegas_fit Objects¶
- class lsqfit.vegas_fit(data, fcn, prior=None, param=None, fit=None, svdcut=1e-12, eps=None, noise=False, **vegasargs)¶
Least-squares fit using Bayesian integrals.
lsqfit.vegas_fit
fits a (nonlinear) functionf(x,p)
(orf(p)
) to datay
using Bayesian integrals over fit parametersp
. Typical usage isvfit = vegas_fit(data=(x,y), fcn=f, prior=prior) print(vfit) print('best-fit parameters =', vfit.p)
The fitter calculates the means and (co)variances of the fit parameters (
vfit.p
) assuming that the parameters are described by a probability density function (PDF) proportional to whereand and . This involves a multi-dimensional integration over the parameter space using
vegas.PDFIntegrator
from thevegas
module (which must be installed separately).vegas
uses adaptive Monte Carlo integration to obtain estimates for the integrals; see its documentation for more information.When the PDF is sufficiently peaked around its maximum, is (usually) well approximated by a quadratic expansion around its minimum, and results obtained from
lsqfit.vegas_fit
will agree with those obtained fromlsqfit.nonlinear_fit
— the latter is the Gaussian approximation to the former. The output fromnonlinear_fit
can often be used to improve significantly the accuracy of the numerical integrals used byvegas_fit
, particularly if the PDF is sharply peaked and there are lots of parameters: for example, by settingparam=fit.p
infit = nonlinear_fit(data=(x,y), fcn=f, prior=prior) vfit = vegas_fit(data=(x,y), fcn=f, prior=prior, param=fit.p) print(vfit)
we direct
vegas_fit
to re-express the integrals overp
in terms of variables that emphasize the region indicated byfit.p
. This facilitates the integration and can greatly reduce the numerical uncertainties in the results. Note that the second line in this code snippet can be written more succinctly asvfit = vegas_fit(fit=fit)
vegas
adapts iteratively to the PDF, averaging results over the iterations. By default it uses 10 iterations to train the integrator on the PDF, and then 10 more, without further adaptation, to estimate the means and covariances of the fit parameters. During the training stage, the integrator remaps the integrals to variables that emphasize regions where the PDF is large, refining the map after each iteration. Integral results from the training stage are often unreliable and so are discarded. (Adaptation is turned off for the latter iterations to provide more robust estimates; see thevegas
documentation for more information.) The number of evaluations off(x,p)
is limited to at most 1000 by default. Both the number of iterations and the number of function evaluations can be specified using parametersnitn
andneval
, respectively: for examplevfit = vegas_fit(fit=fit, nitn=(5,10), neval=10_000)
specifies 5 iterations for adapting to the PDF (training) followed by 10 more for computing
vfit.p
, with at most 10,000 function evaluations per iteration. The number of function evaluations needed depends upon the number of parameters and how sharply peaked the PDF is.Having obtained a fit
vfit
fromvegas_fit
, expectation values with respect to the PDF can be obtained for any functiong(p)
of the parameters, whereg(p)
typically returns an array of numbers or a dictionary whose values are numbers or arrays of numbers, so that multiple expectation values can be computed together:s = vfit.stats(g)
where
s
is an array or dictionary ofgvar.GVar
s giving the means and covariances of the components ofg(p)
with respect to the PDF. Results
agrees well withg(vfit.p)
when the Gaussian approximation is valid (for the PDF andg(p)
), but could be quite different otherwise.vfit.stats
can also calculate other statistical parameters (e.g., skewness) and/or histograms for the distribution ofg(p)
values.- Parameters:
data (dict, array or tuple) –
Data to be fit by
lsqfit.vegas_fit
can have either of the following forms:data = x, y
x
is the independent data that is passed to the fit function with the fit parameters:fcn(x, p)
.y
is a dictionary (or array) ofgvar.GVar
s that encode the means and covariance matrix for the data that is to be fit being fit. The fit function must return a result having the same layout asy
.data = y
y
is a dictionary (or array) ofgvar.GVar
s that encode the means and covariance matrix for the data being fit. There is no independent data so the fit function depends only upon the fit parameters:fit(p)
. The fit function must return a result having the same layout asy
.
Setting
x=False
in the first of these formats implies that the fit function depends only on the fit parameters: that is,fcn(p)
instead offcn(x, p)
. (This is not assumed ifx=None
.) Ignored if parameterfit
is specified.fcn (callable) –
The function to be fit to
data
. It is either a function of the independent datax
and the fit parametersp
(fcn(x, p)
), or a function of just the fit parameters (fcn(p)
) when there is nox
data orx=False
. The function’s return value must have the same layout as they
data (a dictionary or an array). The fit parametersp
are either: 1) a dictionary where eachp[k]
is a single parameter or an array of parameters (any shape); or, 2) a single array of parameters. The layout of the parameters is the same as that of priorprior
if it is specified; otherwise, it is inferred fromparam
.vegas_fit
is usually much faster iffcn
is designed to process a large batch of integration points all at once. See thevegas
documentation onrbatchintegrand
andlbatchintegrand
.prior – A dictionary or array of
gvar.GVar
s representing a priori estimates for all parametersp
used by fit functionfcn(x, p)
(orfcn(p)
). Fit parametersp
are stored in a dictionary or array with the same keys and structure (or shape) asprior
. The default value isNone
;prior
must be defined ifparam
isNone
. Ignored if parameterfit
is specified.param – A dictionary or array of
gvar.GVar
s that specifies the fit parametersp
used byfcn(x,p)
(orfcn(p)
), and indicates where in that parameter space the integrator should focus its attention. Fit parametersp
are stored in a dictionary or array with the same keys and structure (or shape) asparam
(andprior
, if specified).vegas_fit
re-expresses the parameter integrals in terms of variables that emphasize the region of parameter space covered byparam
. Settingparam=None
(the default) is equivalent to settingparam=prior
;prior
must be defined ifparam=None
. Ignored if parameterfit
is specified.fit – Fit results from either
lsqfit.nonlinear_fit
orlsqfit.vegas_fit
. Whenfit
is specified, the data, prior, and fit function are take fromfit
andparam=fit.p
is set. The fit function fromfit
can be replaced by setting thefcn
parameter (for example, to replacefit.fcn
by an equivalent batch function).svdcut (float) – If nonzero, singularities in the correlation matrix for
y
andprior
are regulated usinggvar.regulate()
with an SVD cutoffsvdcut
. This makes the correlation matrices less singular, which can improve the stability and accuracy of a fit. Default issvdcut=1e-12
. Ignored if parameterfit
is specified.eps (float) – If positive, singularities in the correlation matrix for
y
andprior
are regulated usinggvar.regulate()
with cutoffeps
. This makes the correlation matrices less singular, which can improve the stability and accuracy of a fit. Ignored ifsvdcut
is specified (and notNone
). Ignored if parameterfit
is specified.noise (tuple or bool) – If
noise[0]=True
, noise is added to the data means commensurate with the additional uncertainties (if any) introduced by usingsvdcut>0
oreps>0
. Ifnoise[1]=True
, noise is added to the prior means commensurate with the uncertainties in the prior. Noise is useful for testing the quality of a fit (chi2
). Settingnoise=True
is shorthand fornoise=(True, True)
, andnoise=False
meansnoise=(False, False)
(the default). Ignored if parameterfit
is specified.vegasargs – Any additional keyword argments are passed to the integrator,
vegas.PDFIntegrator
. The most important of these arguments are the number ofvegas
interationsnitn
and the maximum numberneval
of integrand evaluations per iteration. Default values for these arenitn=(10,10)
andneval=1000
, where nitn[0] is the number of iterations used to train the integrator to the PDF, and nitn[1] is the number of iterations used to determine the means and covariances of the parameters. (vegas
adapts to the PDF in the first set of (training) iterations; adaptation is turned off for the second.)
Objects of type
lsqfit.nonlinear_fit
have the following attributes:- chi2¶
evaluated at
vfit.p
.fit.chi2 / fit.dof
is usually of order one in good fits. Values much less than one suggest that actual fluctuations in the input data and/or priors might be smaller than suggested by the standard deviations (or covariances) used in the fit.
- dof¶
Number of degrees of freedom in the fit, which equals the number of pieces of data being fit when priors are specified for the fit parameters. Without priors, it is the number of pieces of data minus the number of fit parameters.
- integrator¶
The
vegas.PDFIntegrator
used to do the integrals.
- logBF¶
The logarithm of the probability (density) of obtaining the fit data by randomly sampling the parameter model (priors plus fit function) used in the fit — that is, it is the logarithm of
P(data|model)
. This quantity is useful for comparing fits of the same data to different models, with different priors and/or fit functions. The model with the largest value offit.logBF
is the one preferred by the data. The exponential of the difference infit.logBF
between two models is the ratio of probabilities (Bayes factor) for those models. Differences infit.logBF
smaller than 1 are not very significant.
- p¶
Best-fit parameters from fit in the same format as the prior (array or dictionary containing
gvar.GVar
s). The means and uncertainties (standard deviations/covariances) are calculated from the Bayesian integrals. The uncertainties come from the distribution and from uncertainties invegas
’s estimates of the means (added in quadrature).vfit.p
is output fromvegas.PDFIntegrator.stats()
. It has additional attributes that provide more information about the integrals. See the documentation forvegas.PDFEV
,PDFEVArray
, andvegas.PDFEVDict
for more information.
- pmean¶
An array or dictionary containing the means of the best-fit parameters from fit.
- psdev¶
Standard deviations of the best-fit parameters from fit.
- palt¶
Same as
vfit.p
.
- prior¶
Prior used in the fit. This may differ from the input prior if an SVD cut is set (it is the prior after the SVD cut). It is either a dictionary (
gvar.BufferDict
) or an array (numpy.ndarray
), depending upon the input. EqualsNone
if no prior is specified.
- Q¶
The probability that the
chi**2
from the fit could have been larger, by chance, assuming the best-fit model is correct. Good fits haveQ
values larger than 0.1 or so. Also called the p-value of the fit. The probabilistic intrepretation becomes unreliable if the actual fluctuations in the input data and/or priors are much smaller than suggested by the standard deviations (or covariances) used in the fit (leading to an unusually smallchi**2
).
- residuals¶
An array containing the fit residuals normalized by the corresponding standard deviations. The residuals are projected onto the eigenvectors of the correlation matrix and so should be uncorrelated from each other. The residuals include contributions from both the fit data and the prior. They are related to the the
chi**2
of the fit by:chi2 = sum(fit.residuals**2)
.
- training¶
Results returned by the
vegas.PDFIntegrator
after evaluating the integrals used to train the integrator to the PDF. These results are not included in the final averages.
- correction¶
Sum of all corrections, if any, added to the fit data and prior when
eps>0
orsvdcut>0
.
- svdn¶
Number of eigenmodes of the correlation matrix modified (and/or deleted) when
svdcut>0
.
- time¶
CPU time (in secs) taken by fit.
- x¶
The first field in the input
data
. This is sometimes the independent variable (as in ‘y vs x’ plot), but may be anything. It is set equal toFalse
if thex
field is omitted from the inputdata
. (This also means that the fit function has nox
argument: sof(p)
rather thanf(x,p)
.)
- y¶
Fit data used in the fit, with
gvar.GVar
for each data point. This may differ rom the input data if an SVD cut is used (it isy
after the SVD cut). It is either a dictionary (gvar.BufferDict
) or an array (numpy.ndarray
), depending upon the input.
- nblocks¶
nblocks[s]
equals the number of block-diagonal sub-matrices of they
–prior
covariance matrix that are sizes
-by-s
. This is sometimes useful for debugging.
Additional methods are provided for printing out detailed information about the fit, evaluating
chi**2
:- stats(f, moments=False, histograms=False)¶
Means and standard deviations (and covariances) of function
f(p)
.If
f
is set toNone
or omitted, the means and standard deviations (and covariances) of the fit parameters are recalculated.See documentation for
vegas.PDFIntegrator.stats()
for further options.
- sample(nbatch, mode='rbatch')¶
Generate random samples from PDF used in fit.
See documentation for
vegas.PDFIntegrator.sample()
for more information.- Parameters:
nbatch (int) – The integrator will return at least
nbatch
samples drawn from its PDF. The actual number of samples is the smallest multiple ofself.last_neval
that is equal to or larger thannbatch
. Results are packaged in arrays or dictionaries whose elements have an extra index labeling the different samples in the batch. The batch index is the rightmost index ifmode='rbatch'
; it is the leftmost index ifmode
is'lbatch'
.mode (bool) – Batch mode. Allowed modes are
'rbatch'
or'lbatch'
, corresponding to batch indices that are on the right or the left, respectively. Default ismode='rbatch'
.
- Returns:
A tuple
(wgts,samples)
containing samples drawn from the integrator’s PDF, together with their weightswgts
. The weighted sample points are distributed through parameter space with a density proportional to the PDF.In general,
samples
is either a dictionary or an array depending upon the format oflsqfit.vegas_fit
parameterparam
. For example, ifparam = gv.gvar(dict(s='1.5(1)', v=['3.2(8)', '1.1(4)']))
then
samples['s'][i]
is a sample for parameterp['s']
where indexi=0,1...nbatch(approx)
labels the sample. The corresponding sample forp['v'][d]
, whered=0
or1
, issamples['v'][d, i]
providedmode='rbatch'
, which is the default. (Otherwise it isp['v'][i, d]
, formode='lbatch'
.) The corresponding weight for this sample iswgts[i]
.When
param
is an array,samples
is an array with the same shape plus an extra sample index which is either on the right (mode='rbatch'
, default) or left (mode='lbatch'
).
- format(maxline=0, pstyle='v')¶
Format the output from a
lsqfit.vegas_fit
.See the documentation for
lsqfit.nonlinear_fit.format()
for more information.
- pdf(p)¶
Probability density function for fit parameters
p
.When the
lsqfit.vegas_fit
object is derived from alsqfit.nonlinear_fit
object (vfit = vegas_fit(fit=fit)
), the PDF is . Otherwise it is just . In either case it is not normalized. Divide byvfit.pdfnorm
to normalize it.
- qqplot_residuals(plot=None)¶
Create QQ plot of the fit residuals.
See the documentation for
lsqfit.nonlinear_fit.qqplot_residuals()
for more information.
- plot_residuals(plot=None)¶
Create QQ plot of the fit residuals.
See the documentation for
lsqfit.nonlinear_fit.plot_residuals()
for more information.
Functions¶
- lsqfit.empbayes_fit(z0, fitargs, p0=None, fitter=lsqfit.nonlinear_fit, **minargs)¶
Return
fit
andz
corresponding to the fitlsqfit.nonlinear_fit(**fitargs(z))
that maximizesfit.logGBF
.This function maximizes the logarithm of the Bayes Factor from fit
lsqfit.nonlinear_fit(**fitargs(z))
by varyingz
, starting atz0
. The fit is redone for each value ofz
that is tried, in order to determinefit.logGBF
.The Bayes Factor is proportional to the probability that the data came from the model (fit function and priors) used in the fit.
empbayes_fit()
finds the model or data that maximizes this probability.One application is illustrated by the following code:
import numpy as np import gvar as gv import lsqfit # fit data x = np.array([1., 2., 3., 4.]) y = np.array([3.4422, 1.2929, 0.4798, 0.1725]) # prior prior = gv.gvar(['10(1)', '1.0(1)']) # fit function def fcn(x, p): return p[0] * gv.exp( - p[1] * x) # find optimal dy def fitargs(z): dy = y * z newy = gv.gvar(y, dy) return dict(data=(x, newy), fcn=fcn, prior=prior) fit, z = lsqfit.empbayes_fit(0.1, fitargs) print(fit.format(True))
Here we want to fit data
y
with fit functionfcn
but we don’t know the uncertainties in oury
values. We assume that the relative errors arex
-independent and uncorrelated. We add the errordy
that maximizes the Bayes Factor, as this is the most likely choice. This fit gives the following output:Least Square Fit: chi2/dof [dof] = 0.58 [4] Q = 0.67 logGBF = 7.4834 Parameters: 0 9.44 (18) [ 10.0 (1.0) ] 1 0.9979 (69) [ 1.00 (10) ] Fit: x[k] y[k] f(x[k],p) --------------------------------------- 1 3.442 (54) 3.481 (45) 2 1.293 (20) 1.283 (11) 3 0.4798 (75) 0.4731 (41) 4 0.1725 (27) 0.1744 (23) Settings: svdcut/n = 1e-12/0 tol = (1e-08*,1e-10,1e-10) (itns/time = 3/0.0)
We have, in effect, used the variation in the data relative to the best fit curve to estimate that the uncertainty in each data point is of order 1.6%.
empbayes_fit()
can be used with other fitters: for example, to uselsqfit.vegas_int
instead oflsqfit.nonlinear_fit
(the default) for the fits, replace the next to last line in the code above withfit, z = lsqfit.empbayes_fit(0.1, fitargs, fitter=lsqfit.vegas_fit).
See also
MultiFitter.empbayes_fit()
.- Parameters:
z0 (number, array or dict) – Starting point for search.
fitargs (callable) – Function of
z
that returns a dictionaryargs
containing thelsqfit.nonlinear_fit
arguments corresponding toz
.z
should have the same layout (number, array or dictionary) asz0
.fitargs(z)
can instead return a tuple(args, plausibility)
, whereargs
is again the dictionary forlsqfit.nonlinear_fit
.plausibility
is the logarithm of the a priori probabilitiy thatz
is sensible. Whenplausibility
is provided,lsqfit.empbayes_fit()
maximizes the sumlogGBF + plausibility
. Specifyingplausibility
is a way of steering selections away from completely implausible values forz
.p0 – Fit-parameter starting values for the first fit.
p0
for subsequent fits is set automatically to optimize fitting unless a value is specified byfitargs
.fitter – Fitter to be used. Default is
lsqfit.nonlinear_fit
; also works withlsqfit.vegas_fit
.minargs (dict) – Optional argument dictionary, passed on to
lsqfit.gsl_multiminex
(orlsqfit.scipy_multiminex
), which finds the minimum.
- Returns:
A tuple containing the best fit (object of type
lsqfit.nonlinear_fit
) and the optimal value for parameterz
.
- lsqfit.wavg(datalist, fast=False, prior=None, **fitterargs)¶
Weighted average of
gvar.GVar
s or arrays/dicts ofgvar.GVar
s.The weighted average of
N
gvar.GVar
sxavg = wavg([g1, g2 ... gN])
is what one obtains from a weighted least-squares fit of the collection of
gvar.GVar
s to the one-parameter fit functiondef f(p): return N * [p]
The average is the best-fit value for fit parameter
p
.gvar.GVar
s with smaller standard deviations carry more weight than those with larger standard deviations; and the averages take account of correlations between thegvar.GVar
s.wavg
also works when eachgi
is an array ofgvar.GVar
s or a dictionary whose values aregvar.GVar
s or arrays ofgvar.GVar
s. Corresponding arrays in differentgi
s must have the same dimension, but can have different shapes (the overlapping components are averaged). When thegi
are dictionaries, they need not all have the same keys.Weighted averages can become costly when the number of random samples being averaged is large (100s or more). In such cases it might be useful to set parameter
fast=True
. This causeswavg
to estimate the weighted average by incorporating the random samples one at a time into a running average:result = datalist[0] for di in datalist[1:]: result = wavg([result, di], ...)
This method can be much faster when
len(datalist)
is large, and gives the exact result when there are no correlations between different elements of listdatalist
. The results are approximately correct whendatalist[i]
anddatalist[j]
are correlated fori!=j
.- Parameters:
datalist (list) – The
gvar.GVar
s to be averaged.datalist
is a one-dimensional sequence ofgvar.GVar
s, or of arrays ofgvar.GVar
s, or of dictionaries containinggvar.GVar
s and/or arrays ofgvar.GVar
s. Corresponding arrays in differentdatalist[i]
s must have the same dimension.fast (bool) – If
fast=True
,wavg
averages thedatalist[i]
sequentially. This can be much faster when averaging a large number of sampes but is only approximate if the different elements ofdatalist
are correlated. Default isFalse
.fitterargs (dict) – Additional arguments (e.g.,
svdcut
oreps
) for thelsqfit.nonlinear_fit
fitter used to do the averaging.
- Returns:
The weighted average is returned as a
gvar.GVar
or an array ofgvar.GVar
s or a dictionary ofgvar.GVar
s and arrays ofgvar.GVar
s. Results have the following extra attributes:chi2 -
chi**2
for weighted average.dof - Effective number of degrees of freedom.
- Q - Quality factor Q (or p-value) for fit:
the probability that the
chi**2
could have been larger, by chance, assuming that the data are all Gaussian and consistent with each other. Values smaller than 0.1 or so suggest that the data are not Gaussian or are inconsistent with each other. Also called the p-value.
time - Time required to do average.
- correction - The corrections made to the data
when
svdcut>0
oreps>0
.
fit - Fit returned by
lsqfit.nonlinear_fit
.
- lsqfit.gammaQ(a, x)¶
Return the normalized incomplete gamma function
Q(a,x) = 1-P(a,x)
.Q(a, x) = 1/Gamma(a) * \int_x^\infty dt exp(-t) t ** (a-1) = 1 - P(a, x)
Note that
gammaQ(ndof/2., chi2/2.)
is the probabilty that one could get achi**2
larger thanchi2
withndof
degrees of freedom even if the model used to constructchi2
is correct.
lsqfit.MultiFitter
Classes¶
lsqfit.MultiFitter
provides a framework for building component
systems to fit multiple pieces of data using a set of custom-designed models,
derived from lsqfit.MultiFitterModel
. Each model encapsulates:
a) a particular fit function; b) a recipe for assembling the corresponding fit
data from a dictionary that contains all of the data; and c) a recipe for
assembling a fit prior drawn from a dictionary containing all the priors.
This allows fit problems to be broken down down into more manageable pieces,
which are then aggregated by lsqfit.MultiFitter
into a single fit.
This framework was developed to support the corrfitter
module which is
used to analyze 2-point and 3-point correlators generated in Monte Carlo
simulations of quantum field theories (like QCD). The corrfitter
module provides two models to describe correlators: corrfitter.Corr2
to describe one 2-point correlator, and corrfitter.Corr3
to describe
one 3-point correlator. A typical analysis involves fitting data for a mixture of
2-point and 3-point correlators, with sometimes hundreds of correlators in all.
Each correlator is described by either a Corr2
model or a Corr3
model. A list of models, one for each correlator, is handed
corrfitter.CorrFitter
(derived from lsqfit.MultiFitter
) to
fit the models to the correlator data. The models for different
correlators typically share many fit parameters.
A simpler example of a model is one that encapsulates a linear fit function:
import gvar as gv
import numpy as np
import lsqfit
class Linear(lsqfit.MultiFitterModel):
def __init__(self, datatag, x, intercept, slope):
super(Linear, self).__init__(datatag)
# the independent variable
self.x = np.array(x)
# keys used to find the intercept and slope in a parameter dictionary
self.intercept = intercept
self.slope = slope
def fitfcn(self, p):
try:
return p[self.intercept] + p[self.slope] * self.x
except KeyError:
# slope parameter marginalized/omitted
return len(self.x) * [p[self.intercept]]
def buildprior(self, prior, mopt=None):
" Extract the model's parameters from prior. "
newprior = {}
newprior[self.intercept] = prior[self.intercept]
if mopt is None:
# slope parameter marginalized/omitted if mopt is not None
newprior[self.slope] = prior[self.slope]
return newprior
def builddata(self, data):
" Extract the model's fit data from data. "
return data[self.datatag]
Imagine four sets of data, each corresponding to x=1,2,3,4
, all of which
have the same intercept but different slopes:
data = gv.gvar(dict(
d1=['1.154(10)', '2.107(16)', '3.042(22)', '3.978(29)'],
d2=['0.692(10)', '1.196(16)', '1.657(22)', '2.189(29)'],
d3=['0.107(10)', '0.030(16)', '-0.027(22)', '-0.149(29)'],
d4=['0.002(10)', '-0.197(16)', '-0.382(22)', '-0.627(29)'],
))
To find the common intercept, we define a model for each set of data:
models = [
Linear('d1', x=[1,2,3,4], intercept='a', slope='s1'),
Linear('d2', x=[1,2,3,4], intercept='a', slope='s2'),
Linear('d3', x=[1,2,3,4], intercept='a', slope='s3'),
Linear('d4', x=[1,2,3,4], intercept='a', slope='s4'),
]
This says that data['d3']
, for example, should be fit with function
p['a'] + p['s3'] * np.array([1,2,3,4])
where p
is a dictionary of fit
parameters. Assume that we know a priori that the intercept and slopes are all
order one:
prior = gv.gvar(dict(a='0(1)', s1='0(1)', s2='0(1)', s3='0(1)', s4='0(1)'))
Then we can fit all the data to determine the intercept:
fitter = lsqfit.MultiFitter(models=models)
fit = fitter.lsqfit(data=data, prior=prior)
print(fit)
print('intercept =', fit.p['a'])
The output from this code is:
Least Square Fit:
chi2/dof [dof] = 0.49 [16] Q = 0.95 logGBF = 18.793
Parameters:
a 0.2012 (78) [ 0.0 (1.0) ]
s1 0.9485 (53) [ 0.0 (1.0) ]
s2 0.4927 (53) [ 0.0 (1.0) ]
s3 -0.0847 (53) [ 0.0 (1.0) ]
s4 -0.2001 (53) [ 0.0 (1.0) ]
Settings:
svdcut/n = 1e-12/0 tol = (1e-08*,1e-10,1e-10) (itns/time = 5/0.0)
intercept = 0.2012(78)
Model class Linear
is configured to allow
marginalization of the slope parameter, if desired. Calling
fitter.lsqfit(data=data, prior=prior, mopt=True)
moves the slope
parameters into the data (by subtracting m.x * prior[m.slope]
from the data for each model m
), and does a single-parameter fit for the
intercept:
Least Square Fit:
chi2/dof [dof] = 0.49 [16] Q = 0.95 logGBF = 18.793
Parameters:
a 0.2012 (78) [ 0.0 (1.0) ]
Settings:
svdcut/n = 1e-12/0 tol = (1e-08*,1e-10,1e-10) (itns/time = 4/0.0)
intercept = 0.2012(78)
Marginalization can be useful when fitting large data sets since it reduces the number of fit parameters and simplifies the fit.
Empirical Bayes tuning can be used with a MultiFitter
(see Tuning Priors with the Empirical Bayes Criterion).
Continuing from the example just above, we may be uncertain about the prior for the
intercept. The following code varies the width of that prior to maximize
the Bayes Factor (logGBF
):
def fitargs(z):
prior = gv.gvar(dict(s1='0(1)', s2='0(1)', s3='0(1)', s4='0(1)'))
prior['a'] = gv.gvar(0, np.exp(z)) # np.exp => positive std dev
return dict(prior=prior, data=data, mopt=True)
fit,z = fitter.empbayes_fit(0, fitargs)
print(fit)
print('intercept =', fit.p['a'])
The output shows that the data prefer a prior of 0.0(2)
for the
intercept (not surprisingly):
Least Square Fit:
chi2/dof [dof] = 0.55 [16] Q = 0.92 logGBF = 19.917
Parameters:
a 0.2009 (78) [ 0.00 (20) ]
Settings:
svdcut/n = 1e-12/0 tol = (1e-08*,1e-10,1e-10) (itns/time = 1/0.0)
intercept = 0.2009(78)
The increase in the Bayes Factor, however, is not significant, and the result is almost unchanged. This confirms that the original choice was reasonable.
Another variation is to replace the simultaneous fit of the four models by a chained fit, where one model is fit at a time and its results are fed into the next fit through that fit’s prior. Replacing the fit code by
fitter = lsqfit.MultiFitter(models=models)
fit = fitter.chained_lsqfit(data=data, prior=prior)
# same as fit = fitter.lsqfit(data=data, prior=prior, chained=True)
print(fit.formatall())
print('intercept =', fit.p['a'])
gives the following output:
========== d1
Least Square Fit:
chi2/dof [dof] = 0.32 [4] Q = 0.86 logGBF = 2.0969
Parameters:
a 0.213 (16) [ 0.0 (1.0) ]
s1 0.9432 (82) [ 0.0 (1.0) ]
Settings:
svdcut/n = 1e-12/0 tol = (1e-08*,1e-10,1e-10) (itns/time = 5/0.0)
========== d2
Least Square Fit:
chi2/dof [dof] = 0.58 [4] Q = 0.67 logGBF = 5.3792
Parameters:
a 0.206 (11) [ 0.213 (16) ]
s2 0.4904 (64) [ 0.0 (1.0) ]
Settings:
svdcut/n = 1e-12/0 tol = (1e-08*,1e-10,1e-10) (itns/time = 4/0.0)
========== d3
Least Square Fit:
chi2/dof [dof] = 0.66 [4] Q = 0.62 logGBF = 5.3767
Parameters:
a 0.1995 (90) [ 0.206 (11) ]
s3 -0.0840 (57) [ 0.0 (1.0) ]
Settings:
svdcut/n = 1e-12/0 tol = (1e-08*,1e-10,1e-10) (itns/time = 4/0.0)
========== d4
Least Square Fit:
chi2/dof [dof] = 0.41 [4] Q = 0.81 logGBF = 5.9402
Parameters:
a 0.2012 (78) [ 0.1995 (90) ]
s4 -0.2001 (53) [ 0.0 (1.0) ]
Settings:
svdcut/n = 1e-12/0 tol = (1e-08*,1e-10,1e-10) (itns/time = 4/0.0)
intercept = 0.2012(78)
Note how the value for a
improves with each fit.
Chained fits are most useful
with very large data sets when it is possible to break the data into
smaller, more manageable chunks. There are a variety of options for
organizing the chain of fits; these are discussed in the
MultiFitter.chained_lsqfit()
documentation.
- class lsqfit.MultiFitter(models, mopt=None, ratio=False, fast=True, **fitterargs)¶
Nonlinear least-squares fitter for a collection of models.
Fits collections of data that are modeled by collections of models. Fits can be simultaneous (
lsqfit.MultiFitter.lsqfit()
) or chained (lsqfit.MultiFitter.chained_lsqfit()
).- Parameters:
models – List of models, derived from
lsqfit.MultiFitterModel
, to be fit to the data. Individual models in the list can be replaced by lists of models or tuples of models; see below.mopt (object) – Marginalization options. If not
None
, marginalization is used to reduce the number of fit parameters. Objectmopt
is passed to the models when constructing the prior for a fit; it typically indicates the degree of marginalization (in a model-dependent fashion). Settingmopt=None
implies no marginalization.ratio (bool) – If
True
, implement marginalization using ratios:data_marg = data * fitfcn(prior_marg) / fitfcn(prior)
. IfFalse
(default), implement using differences:data_marg = data + (fitfcn(prior_marg) - fitfcn(prior))
.fast (bool) – Setting
fast=True
(default) strips any variable not required by the fit from the prior. This speeds fits but loses information about correlations between variables in the fit and those that are not. Settingwavg_all=True
can restore some of the correlations, but is somewhat slower.wavg_all (bool) – If
True
andfast=True
, the final result of a chained fit is the weighted average of all the fits in the chain. This can restore correlations lost in the chain becausefast=True
. This step is omitted ifwavg_all=False
orfast=False
. Default isFalse
.fitname (callable or
None
) – Individual fits in a chained fit are assigned default names, constructed from the datatags of the corresponding models, for access and reporting. These names get unwieldy when lots of models are involved. Whenfitname
is notNone
(default), each default namedname
is replaced byfitname(dname)
which should return a string.wavg_kargs (dict) – Keyword arguments for
lsqfit.wavg()
when used to combine results from parallel sub-fits in a chained fit.fitterargs (dict) – Additional arguments for the
lsqfit.nonlinear_fit
object used to do the fits. These can be collected in a dictionary (e.g.,fitterargs=dict(tol=1e-6, maxit=500))
) or listed as separate arguments (e.g.,tol=1e-6, maxit=500
).
- lsqfit(data=None, pdata=None, prior=None, p0=None, chained=False, **kargs)¶
Compute least-squares fit of models to data.
MultiFitter.lsqfit()
fits all of the models together, in a single fit. It returns thelsqfit.nonlinear_fit
object from the fit.To see plots of the fit data divided by the fit function with the best-fit parameters use
fit.show_plots()
This method has optional keyword arguments
save
andview
; see documentation forlsqfit.MultiFitter.show_plots
for more information. Plotting requires modulematplotlib
.To bootstrap a fit, use
fit.bootstrapped_fit_iter(...)
; seelsqfit.nonlinear_fit.bootstrapped_fit_iter()
for more information.- Parameters:
data – Input data. One of
data
orpdata
must be specified but not both.pdata
is obtained fromdata
by collecting the output fromm.builddata(data)
for each modelm
and storing it in a dictionary with keym.datatag
.pdata – Input data that has been processed by the models using
MultiFitter.process_data()
orMultiFitter.process_dataset()
. One ofdata
orpdata
must be specified but not both.prior (dict) – Bayesian prior for fit parameters used by the models.
p0 – Dictionary , indexed by parameter labels, containing initial values for the parameters in the fit. Setting
p0=None
implies that initial values are extracted from the prior. Settingp0="filename"
causes the fitter to look in the file with name"filename"
for initial values and to write out best-fit parameter values after the fit (for the next call toself.lsqfit()
).chained (bool) – Use
MultiFitter.chained_lsqfit()
instead ofMultiFitter.lsqfit()
ifchained=True
. Ignored otherwise. Default ischained=False
.kargs – Arguments that (temporarily) override parameters specified when the
MultiFitter
was created. Can also include additional arguments to be passed through to thelsqfit
fitter.
- chained_lsqfit(data=None, pdata=None, prior=None, p0=None, **kargs)¶
Compute chained least-squares fit of models to data. Equivalent to:
self.lsqfit(data, pdata, prior, p0, chained=True, **kargs).
In a chained fit to models
[s1, s2, ...]
, the models are fit one at a time, with the fit output from one being fed into the prior for the next. This can be much faster than fitting the models together, simultaneously. The final result comes from the last fit in the chain, and includes parameters from all of the models.The most general chain has the structure
[s1, s2, s3 ...]
where eachsn
is one of:A model (derived from
multifitter.MultiFitterModel
).- A tuple
(m1, m2, m3)
of models, to be fit together in a single fit (i.e., simultaneously). Simultaneous fits are useful for closely related models.
- A tuple
- A list
[p1, p2, p3 ...]
where eachpn
is either a model, a tuple of models (see #2), or a dictionary (see #4). The
pn
are fit separately: the fit output from one fit is not fed into the prior of the next (i.e., the fits are effectively in parallel). Results from the separate fits are averaged at the end to provide a single composite result for the collection of fits. Parallel fits are effective (and fast) when the different fits have few or no fit parameters in common.
- A list
- A dictionary that (temporarily) resets default values for
fitter keywords. The new values, specified in the dictionary, apply to subsequent fits in the chain. Any number of such dictionaries can be included in the model chain.
Fit results are returned in a
lsqfit.MultiFitter.chained_nonlinear_fit
objectfit
, which is very similar to anonlinear_fit
object (see documentation for more information). Objectfit
has an extra attributefit.chained_fits
which is an ordered dictionary containing fit results for each link in the chain of fits, indexed by fit names built from the corresponding data tags.To list results from all of the fits in the chain, use
print(fit.formatall())
This method has optional keyword arguments
maxline
,pstyle
, andnline
; see the documentation forlsqfit.nonlinear_fit.format()
for more information.To view plots of each fit use
fit.show_plots()
This method has optional keyword arguments
save
andview
; see documentation forlsqfit.MultiFitter.show_plots
for more information. Plotting requires modulematplotlib
.To bootstrap a fit, use
fit.bootstrapped_fit_iter(...)
; seelsqfit.nonlinear_fit.bootstrapped_fit_iter()
for more information.- Parameters:
data – Input data. One of
data
orpdata
must be specified but not both.pdata
is obtained fromdata
by collecting the output fromm.builddata(data)
for each modelm
and storing it in a dictionary with keym.datatag
.pdata – Input data that has been processed by the models using
MultiFitter.process_data()
orMultiFitter.process_dataset()
. One ofdata
orpdata
must be specified but not both.prior – Bayesian prior for fit parameters used by the models.
p0 –
Dictionary , indexed by parameter labels, containing initial values for the parameters in the fit. Setting
p0=None
implies that initial values are extracted from the prior. Settingp0="filename"
causes the fitter to look in the file with name"filename"
for initial values and to write out best-fit parameter values after the fit (for the next call toself.chained_lsqfit()
). Finally,p0
can be a list containing a differentp0
for each fit in the chain: for example,p0 = [f.pmean for f in fit.chained_fits.values()]
might be a good starting point for the next fit.
kargs – Arguments that override parameters specified when the
MultiFitter
was created. Can also include additional arguments to be passed through to thelsqfit
fitter.
- empbayes_fit(z0, fitargs, p0=None, **minargs)¶
Return fit and
z
corresponding to the fitself.lsqfit(**fitargs(z))
that maximizeslogGBF
.This function maximizes the logarithm of the Bayes Factor from fit
self.lsqfit(**fitargs(z))
by varyingz
, starting atz0
. The fit is redone for each value ofz
that is tried, in order to determinelogGBF
.The Bayes Factor is proportional to the probability that the data came from the model (fit function and priors) used in the fit.
MultiFitter.empbayes_fit()
finds the model or data that maximizes this probability. Seelsqfit.empbayes_fit()
for more information.Include
chained=True
in the dictionary returned byfitargs(z)
if chained fits are desired. See documentation forMultiFitter.lsqfit()
.- Parameters:
z0 (number, array or dict) – Starting point for search.
fitargs (callable) – Function of
z
that returns a dictionaryargs
containing theMultiFitter.lsqfit()
arguments corresponding toz
.z
should have the same layout (number, array or dictionary) asz0
.fitargs(z)
can instead return a tuple(args, plausibility)
, whereargs
is again the dictionary forMultiFitter.lsqfit()
.plausibility
is the logarithm of the a priori probabilitiy thatz
is sensible. Whenplausibility
is provided,MultiFitter.empbayes_fit()
maximizes the sumlogGBF + plausibility
. Specifyingplausibility
is a way of steering selections away from completely implausible values forz
.p0 – Fit-parameter starting values for the first fit.
p0
for subsequent fits is set automatically to optimize fitting unless a value is specified byfitargs
.minargs (dict) – Optional argument dictionary, passed on to
lsqfit.gsl_multiminex
(orlsqfit.scipy_multiminex
), which finds the minimum.
- Returns:
A tuple containing the best fit (a fit object) and the optimal value for parameter
z
.
- set(**kargs)¶
Reset default keyword parameters.
Assigns new default values from dictionary
kargs
to the fitter’s keyword parameters. Keywords for the underlyinglsqfit
fitters can also be included (or grouped together in dictionaryfitterargs
).Returns tuple
(kargs, oldkargs)
wherekargs
is a dictionary containing alllsqfit.MultiFitter
keywords after they have been updated, andoldkargs
contains the original values for these keywords. Usefitter.set(**oldkargs)
to restore the original values.
- static process_data(data, models)¶
Convert
data
to processed data usingmodels
.Data from dictionary
data
is processed by each model in listmodels
, and the results collected into a new dictionarypdata
for use inMultiFitter.lsqfit()
andMultiFitter.chained_lsqft()
.
- static process_dataset(dataset, models, **kargs)¶
Convert
dataset
to processed data usingmodels
.gvar.dataset.Dataset
(or similar dictionary) objectdataset
is processed by each model in listmodels
, and the results collected into a new dictionarypdata
for use inMultiFitter.lsqfit()
andMultiFitter.chained_lsqft()
. Assumes that the models have defined methodMultiFitterModel.builddataset()
. Keyword argumentskargs
are passed on togvar.dataset.avg_data()
when averaging the data.
- static show_plots(fitdata, fitval, x=None, save=False, view='ratio')¶
Show plots comparing
fitdata[k],fitval[k]
for each keyk
infitval
.Assumes
matplotlib
is installed (to make the plots). Plots are shown for one correlator at a time. Press keyn
to see the next correlator; press keyp
to see the previous one; press keyq
to quit the plot and return control to the calling program; press a digit to go directly to one of the first ten plots. Zoom, pan and save using the window controls.There are several different views available for each plot, specified by parameter
view
:view='ratio'
: Data divided by fit (default).view='diff'
: Data minus fit, divided by data’s standard deviation.view='std'
: Data and fit.view='log'
:'std'
with log scale on the vertical axis.view='loglog'
: ‘std’` with log scale on both axes.Press key
v
to cycle through these views; or press keysr
,d
, orl
for the'ratio'
,'diff'
, or'log'
views, respectively.Copies of the plots that are viewed can be saved by setting parameter
save=fmt
wherefmt
is a string used to create file names: the file name for the plot corresponding to keyk
isfmt.format(k)
. It is important that the filename end with a suffix indicating the type of plot file desired: e.g.,fmt='plot-{}.pdf'
.
- static flatten_models(models)¶
Create 1d-array containing all disctinct models from
models
.
lsqfit.MultiFitter
models are derived from the following
class. Methods buildprior
, builddata
, fitfcn
, and
builddataset
are not implemented in this base
class. They need to be overwritten by the derived class (except
for builddataset
which is optional).
- class lsqfit.MultiFitterModel(datatag, ncg=1)¶
Base class for MultiFitter models.
Derived classes must define methods
fitfcn
,buildprior
, andbuilddata
, all of which are described below. In addition they have attributes:- datatag¶
lsqfit.MultiFitter
builds fit data for the correlator by extracting the data labelled bydatatag
(eg, a string) from an input data set (eg, a dictionary). This label is stored in theMultiFitterModel
and must be passed to its constructor. It must be a hashable quantity, like a string or number or tuple of strings and numbers.
- ncg¶
When
ncg>1
, fit data and functions are coarse-grained by breaking them up into bins of ofncg
values and replacing each bin by its average. This can increase the fitting speed, because there is less data, without much loss of precision if the data elements within a bin are highly correlated.
- Parameters:
datatag – Label used to identify model’s data.
ncg (int) – Size of bins for coarse graining (default is
ncg=1
).
- buildprior(prior, mopt=None)¶
Extract fit prior from
prior
.Returns a dictionary containing the part of dictionary
prior
that is relevant to this model’s fit. The code could be as simple as collecting the appropriate pieces: e.g.,def buildprior(self, prior, mopt=None): mprior = gv.BufferDict() model_keys = [...] for k in model_keys: mprior[k] = prior[k] return mprior
where
model_keys
is a list of keys corresponding to the model’s parameters. Supporting non-Gaussian distributions requires a slight modification: e.g.,def buildprior(self, prior, mopt=None): mprior = gv.BufferDict() model_keys = [...] for k in gv.get_dictkeys(prior, model_keys): mprior[k] = prior[k] return mprior
Marginalization involves omitting some of the fit parameters from the model’s prior.
mopt=None
implies no marginalization. Otherwisemopt
will typically contain information about what and how much to marginalize.- Parameters:
prior – Dictionary containing a priori estimates of all fit parameters.
mopt (object) – Marginalization options. Ignore if
None
. Otherwise marginalize fit parameters as specified bymopt
.mopt
can be any type of Python object; it is used only inbuildprior
and is passed through to it unchanged.
- builddata(data)¶
Extract fit data corresponding to this model from data set
data
.The fit data is returned in a 1-dimensional array; the fitfcn must return arrays of the same length.
- Parameters:
data – Data set containing the fit data for all models. This is typically a dictionary, whose keys are the
datatag
s of the models.
- fitfcn(p)¶
Compute fit function fit for parameters
p
.Results are returned in a 1-dimensional array the same length as (and corresponding to) the fit data returned by
self.builddata(data)
.If marginalization is supported,
fitfcn
must work with or without the marginalized parameters.- Parameters:
p – Dictionary of parameter values.
- builddataset(dataset)¶
Extract fit dataset from
gvar.dataset.Dataset
dataset
.The code
import gvar as gv data = gv.dataset.avg_data(m.builddataset(dataset))
that builds data for model
m
should be functionally equivalent toimport gvar as gv data = m.builddata(gv.dataset.avg_data(dataset))
This method is optional. It is used only by
MultiFitter.process_dataset()
.- Parameters:
dataset –
gvar.dataset.Dataset
(or similar dictionary) dataset containing the fit data for all models. This is typically a dictionary, whose keys are thedatatag
s of the models.
References¶
The lsqfit
and gvar
modules were originally created to facilitate statistical analyses of data generated by lattice QCD simulations. Background information about the techniques used in these modules can be found in several articles (on lattice QCD applications):
For a general discussion of Bayesian fitting (and Empirical Bayes) see: G.P. Lepage et al, Nucl.Phys.Proc.Suppl. 106 (2002) 12-20 [hep-lat/0110175].
For a discussion of the underlying analysis in a fit and the meaning of the error budget see Appendix A in: C. Bouchard et al, Phys.Rev. D90 (2014) 054506 [arXiv:1406.2279].
For a discussion of marginalization see the appendix in: C. McNeile et al, Phys.Rev. D82, 034512 (2010) [arXiv:1004.4285]. For another sample application see: K. Hornbostel et al, Phys.Rev. D85 (2012) 031504 [arXiv:1111.1363].
For a discussion of SVD cuts (and goodness-of-fit) see Appendix D in: R.J. Dowdall et al, Phys.Rev. D100 (2019) 9, 094508 [arXiv:1907.01025].
Requirements¶
lsqfit
relies heavily on the gvar
, and numpy
modules.
Also the fitting and minimization routines are from
the Gnu Scientific Library (GSL) and/or the Python scipy
module.