Retrieve Confidence Intervals of fitted TF1

Newbie
1

Hello,

After fitting a TF1 to a dataset, I want to obtain the standard deviation of its predictions.
Following the ConfidenceIntervals.C tutorial and the ROOT::Fit::FitResult documentation, I tried the non-deprecated GetConfidenceIntervals() functions with no success:

  • std::vector<double> GetConfidenceIntervals(double cl, bool norm) returns an empty vector.
  • void GetConfidenceIntervals(unsigned int n, unsigned int stride1, unsigned int stride2, const double * x, double * ci, double cl, bool norm) doesn’t change the value of ci.
  • When trying to re-implement the latter function by myself, I noticed const IModelFunction * FittedFunction() is returning a null, which may be intended behavior according to the following sentence in the class’ documentation:

only with FCN only fits, when fFitFunc=0

The following snippet should reproduce my issue:

#include <iostream>

#include <Fit/BinData.h>
#include <Fit/DataOptions.h>
#include <Fit/Fitter.h>
#include <Fit/FitResult.h>
#include <Fit/Chi2FCN.h>

#include <Math/WrappedMultiTF1.h>

#include <HFitInterface.h>

#include <TF1.h>
#include <TH1D.h>

double Constant(double const *const /*x*/, double const *const par)
{
    return par[0];
}

int Example()
{
    constexpr size_t nevents{10000};
    constexpr size_t nbins{100};
    constexpr double xmin{0}, xmax{15};

    constexpr size_t npar{1};
    constexpr double initial_parameters[npar]{2};

    TF1 *func_gen{new TF1{"func_gen", Constant, xmin, xmax, npar}};
    func_gen->SetParameters(initial_parameters);

    TH1D *hist{new TH1D{"hist", "Random", nbins, xmin, xmax}};
    hist->FillRandom(func_gen->GetName(), nevents);
    hist->Scale(func_gen->Integral(xmin, xmax) / nevents * nbins / (xmax - xmin));

    ROOT::Fit::Fitter fitter{};
    fitter.Config().SetParamsSettings(npar, initial_parameters);

    ROOT::Fit::DataOptions opt{};
    ROOT::Fit::DataRange   range{xmin, xmax};

    ROOT::Fit::BinData bindata{opt, range};
    ROOT::Fit::FillData(bindata, hist);

    TF1                        *func_fit{new TF1{"func_fit", Constant, xmin, xmax, npar}};
    ROOT::Math::WrappedMultiTF1 wf{*func_fit, 1};

    ROOT::Fit::Chi2Function fcn{bindata, wf};

    fitter.FitFCN(npar, fcn, 0, bindata.Size(), true);
    const auto result{fitter.Result()};

    func_fit->SetFitResult(result);

    func_gen->SetLineColor(kRed);
    hist->GetListOfFunctions()->Add(func_gen);
    func_fit->SetLineColor(kBlue);
    hist->GetListOfFunctions()->Add(func_fit);
    hist->Draw();

    result.Print(std::cout);

    std::cout << "result.GetConfidenceIntervals().size() = " << result.GetConfidenceIntervals().size() << std::endl;

    double x{1.}, err{42.};
    result.GetConfidenceIntervals(1, 1, 1, &x, &err);
    std::cout << "err [result.GetConfidenceIntervals(1, 1, 1, &x, &err)] = " << err << std::endl;

    std::cout << "result.FittedFunction() = " << result.FittedFunction() << std::endl;

    return 0;
}

In this case, replacing the use of ROOT::Fit with hist->Fit(func_fit) and (TVirtualFitter::GetFitter())->GetConfidenceIntervals(grint) is likely enough to get the results I want.
However my particular application is is derived from the combinedFit.C tutorial for fitting two functions with common parameters.

Best regards,
Bruno Oliveira

ROOT Version: 6.26/06 (pacman -S root)
Platform: Arch Linux
Compiler: Cling

2

Hi @BrunoOliveira,

welcome back to the root forum. Maybe @moneta could take a look?

Cheers,
Marta

3

Let’s ping @moneta again, just in case.

Cheers,
J.

4

Hi,

If you are using the fitter by providing yourself t he FCN z(chi2 function) it will not work the current functionality we have to get the confidence levels, since the Fitter/FitResult does not know on which function should compute them.
You can probably compute them yourself using the code in ROOT: math/mathcore/src/FitResult.cxx Source File using the correct fitted function.

Lorenzo

5

Hi Lorenzo,

Thank you for your insight.

From my understanding, the correct fitted function in the context of the above example is:
ROOT::Math::WrappedMultiTF1 wf{*func_fit, 1};
and therefore the proper implementation would be:

#include <iostream>

#include <Fit/BinData.h>
#include <Fit/DataOptions.h>
#include <Fit/Fitter.h>
#include <Fit/FitResult.h>
#include <Fit/Chi2FCN.h>

#include <Math/OneDimFunctionAdapter.h>
#include <Math/QuantFunc.h>
#include <Math/RichardsonDerivator.h>
#include <Math/WrappedMultiTF1.h>

#include <HFitInterface.h>

#include <TF1.h>
#include <TH1D.h>

double GetConfidenceInterval(const ROOT::Fit::FitResult &result, const ROOT::Math::IParamMultiFunction &wf,
                             const double x, const double cl = 0.683)
{
    const unsigned int ndf{result.Ndf()};
    const double       chi2{result.Chi2()};

    double     corrFactor{1.};
    const bool norm{result.NormalizedErrors() && !(chi2 <= 0 || ndf == 0)};
    if (norm)
        corrFactor = TMath::StudentQuantile(0.5 + cl / 2, ndf) * std::sqrt(chi2 / ndf);
    else
        corrFactor = ROOT::Math::normal_quantile(0.5 + cl / 2, 1);

    unsigned int        npar{wf.NPar()};
    std::vector<double> grad(npar);
    std::vector<double> vsum(npar);

    ROOT::Math::RichardsonDerivator d;
    for (unsigned int ipar{0}; ipar < npar; ++ipar) {
        if (result.IsParameterFixed(ipar)) {
            grad[ipar] = 0;
        } else {
            ROOT::Math::OneDimParamFunctionAdapter<const ROOT::Math::IParamMultiFunction &> fadapter(
                wf, &x, &result.Parameters().front(), ipar);
            d.SetFunction(fadapter);

            if (result.Error(ipar) > 0)
                d.SetStepSize(std::max(result.Error(ipar) * 1.E-5, 1.E-15));
            else
                d.SetStepSize(std::min(std::max(result.Parameter(ipar) * 1.E-5, 1.E-15), 0.0001));

            grad[ipar] = d(result.Parameter(ipar));
        }
    }

    vsum.assign(npar, 0.0);
    for (unsigned int ipar{0}; ipar < npar; ++ipar) {
        for (unsigned int jpar{0}; jpar < npar; ++jpar) {
            vsum[ipar] += result.CovMatrix(ipar, jpar) * grad[jpar];
        }
    }

    double r2{0};
    for (unsigned int ipar{0}; ipar < npar; ++ipar) {
        r2 += grad[ipar] * vsum[ipar];
    }

    return corrFactor * std::sqrt(r2);
}

double Constant(double const *const /*x*/, double const *const par)
{
    return par[0];
}

int Example()
{
    constexpr size_t nevents{10000};
    constexpr size_t nbins{100};
    constexpr double xmin{0}, xmax{15};

    constexpr size_t npar{1};
    constexpr double initial_parameters[npar]{2};

    TF1 *func_gen{new TF1{"func_gen", Constant, xmin, xmax, npar}};
    func_gen->SetParameters(initial_parameters);

    TH1D *hist{new TH1D{"hist", "Random", nbins, xmin, xmax}};
    hist->FillRandom(func_gen->GetName(), nevents);
    hist->Scale(func_gen->Integral(xmin, xmax) / nevents * nbins / (xmax - xmin));

    ROOT::Fit::Fitter fitter{};
    fitter.Config().SetParamsSettings(npar, initial_parameters);

    ROOT::Fit::DataOptions opt{};
    ROOT::Fit::DataRange   range{xmin, xmax};

    ROOT::Fit::BinData bindata{opt, range};
    ROOT::Fit::FillData(bindata, hist);

    TF1                        *func_fit{new TF1{"func_fit", Constant, xmin, xmax, npar}};
    ROOT::Math::WrappedMultiTF1 wf{*func_fit, 1};

    ROOT::Fit::Chi2Function fcn{bindata, wf};

    fitter.FitFCN(npar, fcn, 0, bindata.Size(), true);
    const auto result{fitter.Result()};

    func_fit->SetFitResult(result);

    func_gen->SetLineColor(kRed);
    hist->GetListOfFunctions()->Add(func_gen);
    func_fit->SetLineColor(kBlue);
    hist->GetListOfFunctions()->Add(func_fit);
    hist->Draw();

    result.Print(std::cout);

    std::cout << "result.GetConfidenceIntervals().size() = " << result.GetConfidenceIntervals().size() << std::endl;

    double x{1.}, err{42.};
    result.GetConfidenceIntervals(1, 1, 1, &x, &err);
    std::cout << "err [result.GetConfidenceIntervals(1, 1, 1, &x, &err)] = " << err << std::endl;

    std::cout << "result.FittedFunction() = " << result.FittedFunction() << std::endl;

    std::cout << "GetConfidenceInterval(result, wf, x) = " << GetConfidenceInterval(result, wf, x) << std::endl;

    return 0;
}

Is this correct, or am I missing something?

Bruno