dpfit

Least-squares fitting utilities developed by A. Merand.

bootstrap(func, x, params, y, err=None, fitOnly=None, verbose=False, doNotFit=None, epsfcn=1e-07, ftol=1e-05, fullOutput=True, normalizedUncer=True, follow=None, Nboot=None)

Bootstrap least-squares fits by resampling data indices.

Parameters:

Name	Type	Description	Default
func	callable	Model called as func(x, params).	required
x	array_like	Independent-variable values to resample.	required
params	dict[str, float]	Initial parameter values.	required
y	array_like	Observed values to resample with x.	required
err	array_like or ndarray	Errors or a two-dimensional error covariance matrix. This is passed unchanged to every resampled fit.	None
fitOnly	list of str	Parameters to fit. All parameters are fitted when omitted.	None
verbose	bool or int	Whether to print fit information. Bootstrap fits are always silent.	False
doNotFit	list of str	Parameters to exclude from fitting.	None
epsfcn	float	Step length for the least-squares solver.	1e-7
ftol	float	Relative-error tolerance for the least-squares solver.	1e-5
fullOutput	bool	Whether to return fit metadata. Bootstrap fits always return full output.	True
normalizedUncer	bool	Whether to scale uncertainties by reduced chi-squared. Bootstrap fits always use normalized uncertainties.	True
follow	list of str	Parameters to report while fitting. Bootstrap fits do not report progress.	None
Nboot	int	Number of resampled fits. Defaults to ten times len(x).	None

Returns:

Type	Description
list of dict	Full-output dictionaries for the original fit followed by Nboot fits based on independently resampled indices.

Source code in src/amical/dpfit.py

def bootstrap(
    func,
    x,
    params,
    y,
    err=None,
    fitOnly=None,
    verbose=False,
    doNotFit=None,
    epsfcn=1e-7,
    ftol=1e-5,
    fullOutput=True,
    normalizedUncer=True,
    follow=None,
    Nboot=None,
):
    """Bootstrap least-squares fits by resampling data indices.

    Parameters
    ----------
    func : callable
        Model called as ``func(x, params)``.
    x : array_like
        Independent-variable values to resample.
    params : dict[str, float]
        Initial parameter values.
    y : array_like
        Observed values to resample with ``x``.
    err : array_like or numpy.ndarray, optional
        Errors or a two-dimensional error covariance matrix. This is passed
        unchanged to every resampled fit.
    fitOnly : list of str, optional
        Parameters to fit. All parameters are fitted when omitted.
    verbose : bool or int, default=False
        Whether to print fit information. Bootstrap fits are always silent.
    doNotFit : list of str, optional
        Parameters to exclude from fitting.
    epsfcn : float, default=1e-7
        Step length for the least-squares solver.
    ftol : float, default=1e-5
        Relative-error tolerance for the least-squares solver.
    fullOutput : bool, default=True
        Whether to return fit metadata. Bootstrap fits always return full
        output.
    normalizedUncer : bool, default=True
        Whether to scale uncertainties by reduced chi-squared. Bootstrap fits
        always use normalized uncertainties.
    follow : list of str, optional
        Parameters to report while fitting. Bootstrap fits do not report
        progress.
    Nboot : int, optional
        Number of resampled fits. Defaults to ten times ``len(x)``.

    Returns
    -------
    list of dict
        Full-output dictionaries for the original fit followed by ``Nboot``
        fits based on independently resampled indices.
    """
    if doNotFit is None:
        doNotFit = []

    if Nboot is None:
        Nboot = 10 * len(x)
    # first fit is the "normal" one
    fits = [
        leastsqFit(
            func,
            x,
            params,
            y,
            err=err,
            fitOnly=fitOnly,
            verbose=False,
            doNotFit=doNotFit,
            epsfcn=epsfcn,
            ftol=ftol,
            fullOutput=True,
            normalizedUncer=True,
        )
    ]
    for _ in range(Nboot):
        s = np.int_(len(x) * np.random.rand(len(x)))
        fits.append(
            leastsqFit(
                func,
                x[s],
                params,
                y[s],
                err=err,
                fitOnly=fitOnly,
                verbose=False,
                doNotFit=doNotFit,
                epsfcn=epsfcn,
                ftol=ftol,
                fullOutput=True,
                normalizedUncer=True,
            )
        )
    return fits

leastsqFit(func, x, params, y, err=None, fitOnly=None, verbose=False, doNotFit=None, epsfcn=1e-07, ftol=1e-05, fullOutput=True, normalizedUncer=True, follow=None)

Fit a callable to data with least squares.

Parameters:

Name	Type	Description	Default
func	callable	Model called as func(x, params).	required
x		Independent and observed values.	required
params	dict	Initial parameter values.	required
y		Independent and observed values.	required
err	array - like or ndarray	Errors or a two-dimensional error covariance matrix.	None
fitOnly	list of str	Parameters to fit; all parameters are fitted by default.	None
verbose	bool or int	Whether to print fit information.	False
doNotFit	list of str	Parameters to exclude from fitting.	None
epsfcn	float	Least-squares solver tolerances.	1e-7, 1e-5
ftol	float	Least-squares solver tolerances.	1e-7, 1e-5
fullOutput	bool	Whether to return fit metadata.	True
normalizedUncer	bool	Whether to scale uncertainties by reduced chi-squared.	True
follow	list of str	Parameters to report while fitting.	None

Returns:

Type	Description
dict	Best-fit parameters or, when fullOutput is true, fit parameters, uncertainties, reduced chi-squared, model, covariance, and metadata.

Source code in src/amical/dpfit.py

def leastsqFit(
    func,
    x,
    params,
    y,
    err=None,
    fitOnly=None,
    verbose=False,
    doNotFit=None,
    epsfcn=1e-7,
    ftol=1e-5,
    fullOutput=True,
    normalizedUncer=True,
    follow=None,
):
    """Fit a callable to data with least squares.

    Parameters
    ----------
    func : callable
        Model called as func(x, params).
    x: array-like
        Independent and observed values.
    params : dict
        Initial parameter values.
    y: array-like
        Independent and observed values.
    err : array-like or numpy.ndarray, optional
        Errors or a two-dimensional error covariance matrix.
    fitOnly : list of str, optional
        Parameters to fit; all parameters are fitted by default.
    verbose : bool or int, default=False
        Whether to print fit information.
    doNotFit : list of str, optional
        Parameters to exclude from fitting.
    epsfcn, ftol : float, default=1e-7, 1e-5
        Least-squares solver tolerances.
    fullOutput : bool, default=True
        Whether to return fit metadata.
    normalizedUncer : bool, default=True
        Whether to scale uncertainties by reduced chi-squared.
    follow : list of str, optional
        Parameters to report while fitting.

    Returns
    -------
    dict
        Best-fit parameters or, when fullOutput is true, fit parameters,
        uncertainties, reduced chi-squared, model, covariance, and metadata.
    """
    import scipy.optimize

    if doNotFit is None:
        doNotFit = []

    # -- fit all parameters by default
    if fitOnly is None:
        if len(doNotFit) > 0:
            fitOnly = [x for x in list(params.keys()) if x not in doNotFit]
        else:
            fitOnly = list(params.keys())
        fitOnly.sort()  # makes some display nicer

    # -- build fitted parameters vector:
    pfit = [params[k] for k in fitOnly]

    # -- built fixed parameters dict:
    pfix = {}
    for k in list(params.keys()):
        if k not in fitOnly:
            pfix[k] = params[k]
    if verbose:
        print(f"[dpfit] {len(fitOnly)} FITTED parameters:", fitOnly)
    # -- actual fit
    plsq, cov, info, mesg, ier = scipy.optimize.leastsq(
        _fitFunc,
        pfit,
        args=(
            fitOnly,
            x,
            y,
            err,
            func,
            pfix,
            verbose,
            follow,
        ),
        full_output=True,
        epsfcn=epsfcn,
        ftol=ftol,
    )
    if isinstance(err, np.ndarray) and len(err.shape) == 2:
        print(cov)

    # -- best fit -> agregate to pfix
    for i, k in enumerate(fitOnly):
        pfix[k] = plsq[i]

    # -- reduced chi2
    model = func(x, pfix)
    tmp = _fitFunc(plsq, fitOnly, x, y, err, func, pfix)

    try:
        chi2 = (np.array(tmp) ** 2).sum()
    except Exception:
        chi2 = 0.0
        for x in tmp:
            chi2 += np.sum(x**2)
    reducedChi2 = chi2 / float(
        np.sum([1 if np.isscalar(i) else len(i) for i in tmp]) - len(pfit) + 1
    )
    # print(chi2, reducedChi2, float(np.sum([1 if np.isscalar(i) else
    #                                 len(i) for i in tmp])-len(pfit)+1))
    if not np.isscalar(reducedChi2):
        reducedChi2 = np.mean(reducedChi2)

    # -- uncertainties:
    uncer = {}
    for k in list(pfix.keys()):
        if k not in fitOnly:
            uncer[k] = 0  # not fitted, uncertatinties to 0
        else:
            i = fitOnly.index(k)
            if cov is None:
                uncer[k] = -1
            else:
                uncer[k] = np.sqrt(np.abs(np.diag(cov)[i]))
                if normalizedUncer:
                    uncer[k] *= np.sqrt(reducedChi2)

    if verbose:
        print("-" * 30)
        print("REDUCED CHI2=", reducedChi2)
        print("-" * 30)
        if normalizedUncer:
            print("(uncertainty normalized to data dispersion)")
        else:
            print("(uncertainty assuming error bars are correct)")
        tmp = list(pfix.keys())
        tmp.sort()
        maxLength = np.max(np.array([len(k) for k in tmp]))
        format_ = "'%s':"
        # -- write each parameter and its best fit, as well as error
        # -- writes directly a dictionnary
        print("")  # leave some space to the eye
        for ik, k in enumerate(tmp):
            padding = " " * (maxLength - len(k))
            formatS = format_ + padding
            if ik == 0:
                formatS = "{" + formatS
            if uncer[k] > 0:
                ndigit = -int(np.log10(uncer[k])) + 3
                print(formatS % k, round(pfix[k], ndigit), ",", end=" ")
                print("# +/-", round(uncer[k], ndigit))
            elif uncer[k] == 0:
                if isinstance(pfix[k], str):
                    print(formatS % k, "'" + pfix[k] + "'", ",")
                else:
                    print(formatS % k, pfix[k], ",")
            else:
                print(formatS % k, pfix[k], ",", end=" ")
                print("# +/-", uncer[k])
        print("}")  # end of the dictionnary
        try:
            if verbose > 1:
                print("-" * 3, "correlations:", "-" * 15)
                N = np.max([len(k) for k in fitOnly])
                N = min(N, 20)
                N = max(N, 5)
                sf = "%" + str(N) + "s"
                print(" " * N, end=" ")
                for k2 in fitOnly:
                    print(sf % k2, end=" ")
                print("")
                sf = "%-" + str(N) + "s"
                for k1 in fitOnly:
                    i1 = fitOnly.index(k1)
                    print(sf % k1, end=" ")
                    for k2 in fitOnly:
                        i2 = fitOnly.index(k2)
                        if k1 != k2:
                            print(
                                ("%" + str(N) + ".2f")
                                % (cov[i1, i2] / np.sqrt(cov[i1, i1] * cov[i2, i2])),
                                end=" ",
                            )
                        else:
                            print(" " * (N - 4) + "-" * 4, end=" ")
                    print("")
                print("-" * 30)
        except Exception:
            pass
    # -- result:
    if fullOutput:
        if normalizedUncer:
            try:
                cov *= reducedChi2
            except Exception:
                pass
        try:
            cor = np.sqrt(np.diag(cov))
            cor = cor[:, None] * cor[None, :]
            cor = cov / cor
        except Exception:
            cor = None

        pfix = {
            "best": pfix,
            "uncer": uncer,
            "chi2": reducedChi2,
            "model": model,
            "cov": cov,
            "fitOnly": fitOnly,
            "info": info,
            "cor": cor,
        }
    return pfix