Distribution from createProjection has erroneus plot

Ako_b · September 10, 2025, 6:25pm

Hello!
I am learning RooFit and have encountered a problem when doing a sanity check with this example:

Following script produces plot of compound pdf (exponential distribution with rate following gamma distribution) in blue and Lomax pdf in red. They should coincide.

void roofit_compound_model_expo() {
    // Observable
    RooRealVar x("x", "x", 1, 3);

    RooRealVar lambda("lambda", "rate", 1, 0, 100);

    RooGenericPdf expo("expo", "expo",
                       "lambda * exp(-lambda * x)",
                       RooArgList(x, lambda));

    // Gamma parameters
    RooRealVar alpha("alpha", "shape", 3.0, 0.1, 10);
    RooRealVar theta("theta", "scale", 1.0, 0.1, 10);

    RooGenericPdf gammaDistribution("gammaDistribution", "gammaDistribution",
                                    "pow(lambda/theta, alpha - 1) * exp(-lambda/theta) / theta / TMath::Gamma(alpha)",
                                    RooArgList(lambda, alpha, theta));


    // Compound model: product of exponential conditional on rate and gamma
    RooProdPdf compoundModel("compoundModel", "expo with gamma rate",
                             RooArgList(expo, gammaDistribution));

    // Create marginal distribution by integrating out lambda
    auto marginal_x = compoundModel.createProjection(RooArgSet(lambda));

    // Lomax distribution as analytical model for comparison
    // theta = 1 / lambda_wiki, where lambda_wiki is parameter from Wiki article 
    RooGenericPdf lomax("lomax", "lomax",
                        "alpha * theta * pow(1 + x * theta, -alpha-1)", RooArgList(x, alpha, theta));

    RooPlot *xFrame = x.frame();
    marginal_x->plotOn(xFrame);
    lomax.plotOn(xFrame, RooFit::LineColor(kRed));

    xFrame->Draw();
}

However, I observe some artifact.

I don’t understand why compound pdf goes to zero abruptly. How to fix this?

Setup

ROOT v6.36.04
Built for linuxx8664gcc on Sep 04 2025, 09:14:25
From tags/6-36-04@6-36-04
With clang version 19.1.7 (AlmaLinux OS Foundation 19.1.7-2.el9)
Binary directory: /opt/root-6.36.04-debug/bin

Wile_E_Coyote · September 11, 2025, 7:32pm

To me, it looks like a usual/known ROOT integrator’s misbehaviour.

{
  RooRealVar x("x", "x", 1, 3);
  RooRealVar lambda("lambda", "rate", 1, 0, 100);
  RooGenericPdf expo("expo", "expo",
                     "lambda * exp(-lambda * x)",
                     RooArgList(x, lambda));
  auto marginal_x = expo.createProjection(RooArgSet(lambda));
  RooPlot *xFrame = x.frame();
  marginal_x->plotOn(xFrame);
  xFrame->Draw();
}

StephanH · September 15, 2025, 3:11pm

@jonas, could you have a look here?

Ako_b · September 15, 2025, 7:38pm

I’ve investigated the issue for a while. TL;DR: RooRombergIntegrator is used for integration, and I suspect that it is buggy.

Test setup

MyExpo.cxx (463 Bytes)
MyExpo.hxx (1.1 KB)
MyLinkDef.hxx (161 Bytes)
roofit_test_create_projection.cxx (911 Bytes)

ROOT v6.36.04
Built for linuxx8664gcc on Sep 04 2025, 09:14:25
From tags/6-36-04@6-36-04
With clang version 19.1.7 (AlmaLinux OS Foundation 19.1.7-2.el9)
Binary directory: /opt/root-6.36.04-debug/bin

Compilation

rootcling -f G_MyExpo.cxx -c MyExpo.hxx MyLinkDef.hxx
g++ -std=c++17 -g -O0 -fPIE -Wall -Wextra -pedantic -Werror \
        -o roofit_test_create_projection \
        roofit_test_create_projection.cxx MyExpo.cxx G_MyExpo.cxx \
        $(root-config --libs --cflags) -lRooFit -lRooFitCore

Result

$ ./roofit_test_create_projection
[#1] INFO:NumericIntegration -- RooRealIntegral::init(expo_Int[lambda]) using numeric integrator RooIntegrator1D to calculate Int(lambda)
---
[#1] INFO:NumericIntegration -- RooRealIntegral::init(expo_Int[lambda]_Norm[lambda]) using numeric integrator RooIntegrator1D to calculate Int(lambda)
         x = 2.33755
marginal_x = 1
---
         x = 2.33765
marginal_x = 7150.76
---

Debug attempt

gdb reveals RooRombergIntegrator under the hood

$ gdb roofit_test_create_projection
(gdb) start
(gdb) break RooRombergIntegrator.cxx:163
(gdb) continue
(gdb) backtrace
#0  RooFit::Detail::integrate1d(std::function<double (double)>, bool, int, int, int, double, double, bool, double, double, std::__ROOT::span<double>, std::__ROOT::span<double>) (func=..., doTrapezoid=true, maxSteps=20, minStepsZero=999, 
    fixSteps=0, epsAbs=9.9999999999999995e-08, epsRel=9.9999999999999995e-08, doExtrap=true, xmin=0, xmax=100, hArr=..., sArr=...)
    at /opt/root-6.36.04-source/roofit/roofitcore/src/RooRombergIntegrator.cxx:167
#1  0x00007ffff41b407b in RooRombergIntegrator::integral (this=0x24a7620, iDim=0, nSeg=1, wksp=...) at /opt/root-6.36.04-source/roofit/roofitcore/src/RooRombergIntegrator.cxx:540
#2  0x00007ffff41b3c30 in RooRombergIntegrator::integral (this=0x24a7620, yvec=0x0) at /opt/root-6.36.04-source/roofit/roofitcore/src/RooRombergIntegrator.cxx:488
#3  0x00007ffff3f511a5 in RooAbsIntegrator::calculate (this=0x24a7620, yvec=0x0) at /opt/root-6.36.04-source/roofit/roofitcore/src/RooAbsIntegrator.cxx:54
#4  0x00007ffff419a102 in RooRealIntegral::integrate (this=0x24a61c0) at /opt/root-6.36.04-source/roofit/roofitcore/src/RooRealIntegral.cxx:977
#5  0x00007ffff419a012 in RooRealIntegral::sum (this=0x24a61c0) at /opt/root-6.36.04-source/roofit/roofitcore/src/RooRealIntegral.cxx:963
#6  0x00007ffff4198f05 in RooRealIntegral::evaluate (this=0x24a61c0) at /opt/root-6.36.04-source/roofit/roofitcore/src/RooRealIntegral.cxx:846
#7  0x00007ffff3f6c6b1 in RooAbsReal::traceEval (this=0x24a61c0) at /opt/root-6.36.04-source/roofit/roofitcore/src/RooAbsReal.cxx:320
#8  0x00007ffff4198b75 in RooRealIntegral::getValV (this=0x24a61c0, nset=0x0) at /opt/root-6.36.04-source/roofit/roofitcore/src/RooRealIntegral.cxx:799
#9  0x000000000040bc31 in RooAbsReal::getVal (this=0x24a61c0, normalisationSet=0x0) at /opt/root-6.36.04-debug/include/RooAbsReal.h:117
#10 0x00007ffff3f5ad08 in RooAbsPdf::syncNormalization (this=0x7fffffffc550, nset=0x24441a0, adjustProxies=true) at /opt/root-6.36.04-source/roofit/roofitcore/src/RooAbsPdf.cxx:541
#11 0x00007ffff3f59a64 in RooAbsPdf::getValV (this=0x7fffffffc550, nset=0x24441a0) at /opt/root-6.36.04-source/roofit/roofitcore/src/RooAbsPdf.cxx:334
#12 0x000000000040bc31 in RooAbsReal::getVal (this=0x7fffffffc550, normalisationSet=0x24441a0) at /opt/root-6.36.04-debug/include/RooAbsReal.h:117
#13 0x00007ffff4195656 in RooRealIntegral::RooRealIntegral (this=0x24354d0, name=0x24354a0 "expo_Int[lambda]_Norm[lambda]", title=0x241d370 "Integral of expo", function=..., depList=..., funcNormSet=0x7fffffffc280, config=0x0, 
    rangeName=0x0) at /opt/root-6.36.04-source/roofit/roofitcore/src/RooRealIntegral.cxx:519
#14 0x00007ffff3f8eaa1 in std::make_unique<RooRealIntegral, char const*, char const*, RooAbsReal const&, RooArgSet&, RooArgSet const*&, RooNumIntConfig const*&, char const*&> (__args=@0x7fffffffbb70: 0x0, __args=@0x7fffffffbb70: 0x0, 
    __args=@0x7fffffffbb70: 0x0, __args=@0x7fffffffbb70: 0x0, __args=@0x7fffffffbb70: 0x0, __args=@0x7fffffffbb70: 0x0, __args=@0x7fffffffbb70: 0x0)
    at /opt/rh/gcc-toolset-14/root/usr/lib/gcc/x86_64-redhat-linux/14/../../../../include/c++/14/bits/unique_ptr.h:1077
#15 0x00007ffff3f6eaa9 in RooAbsReal::createIntObj (this=0x7fffffffc550, iset2=..., nset2=0x7fffffffc280, cfg=0x0, rangeName=0x0) at /opt/root-6.36.04-source/roofit/roofitcore/src/RooAbsReal.cxx:634
#16 0x00007ffff3f6dcb1 in RooAbsReal::createIntegral (this=0x7fffffffc550, iset=..., nset=0x7fffffffc280, cfg=0x0, rangeName=0x0) at /opt/root-6.36.04-source/roofit/roofitcore/src/RooAbsReal.cxx:557
#17 0x00007ffff41859d8 in RooProjectedPdf::getProjection (this=0x2430c30, iset=0x7fffffffd220, nset=0x0, rangeName=0x0, code=@0x7fffffffc3d8: -137947173)
    at /opt/root-6.36.04-source/roofit/roofitcore/src/RooProjectedPdf.cxx:126
#18 0x00007ffff4185647 in RooProjectedPdf::RooProjectedPdf (this=0x2430c30, name=0x241c430 "expo_Proj[lambda]", title=0x241c430 "expo_Proj[lambda]", _intpdf=..., intObs=...)
    at /opt/root-6.36.04-source/roofit/roofitcore/src/RooProjectedPdf.cxx:64
#19 0x00007ffff3f63ca8 in RooAbsPdf::createProjection (this=0x7fffffffc550, iset=...) at /opt/root-6.36.04-source/roofit/roofitcore/src/RooAbsPdf.cxx:2450
#20 0x000000000040b6ae in main () at roofit_test_create_projection.cxx:13

Discussion

Apparently RooRombergIntegrator implementation, which is used by default, has either numerical stability issues or some bugs. I’ve narrowed down the problem to this class by inserting its code into another test setup, which I can send if needed.

RooRombergIntegrator has integrate1d as main routine, which iteratively calculates integral sum via addTrapezoids for finer and finer grids. It also uses extrapolate to improve accuracy using previously calculated results.

After a glance at source code and documentation I’ve noted several weird or unclear things:

Reference page for RooRombergIntegrator is available only for 6.34 as of now, but the class is still present in 6.36.
The class implementation uses std::span, which was introduced in C++20, however, I’ve built my debug ROOT version with C++17 without problems.
extrapolate function looks most suspicious to me, as I don’t understand its code completely, and I didn’t notice anything wrong in other functions. Documentation suggests that Richardson extrapolation is applied only fixed amount of times, contrary to Wiki implementation which applies it more and more times as grid step shrinks. In theory this fixed number of extrapolations reduces accuracy but might avoid stability issues with higher order extrapolation.

I hope that this information will be useful!

Ako_b · September 17, 2025, 9:04am

I think I’ve figured out the problem. The implementation of RooRombergIntegrator seems to be correct. It just stops prematurely after the first extrapolation step.

I’ve plotted how approximations of \int_0^{100} x e^{-2.33875 x} dx with trapezoidal rule behave with smaller and smaller grid size.

Values on the right tail correspond to coarse grid which doesn’t sample well the region \approx (0, 6) where the integrand’s value is concentrated. Conceptually, RooRombergIntegrator extrapolates this plot to the left starting from the five points on the right, and arrives at very small value with seemingly good error estimate.

It doesn’t mean extrapolation procedure is faulty, as it must be robust if trapezoidal sum behaves as a convex or concave function of grid step, which is true for sufficiently small step values. It is not surprising that it fails to capture this S-like shape using only information from its tail.

By the way, it seems that trapezoidal sum is a logarithmically convex function of grid step. Extrapolating logarithms of sums might be more robust in case with coarse grids and near zero functions.

I guess that tuning RooRombergIntegrator or switching to other integrator might fix the initial problem as is.
Otherwise, I see several solutions

Make class with analyticalIntegral implementation
Use non-rectangular domains (in the first post, different integration limits for \lambda depending on value of x)
Use Monte Carlo integration (source: RooFit Users Manual v2.91, Chapter 6, Section Conditional probability density functions)

StephanH · September 17, 2025, 10:19am

Hello @Ako_b,

nice analysis of the integrator! I can see that it can stop too early for certain “misbehaving” functions (in the sense that the series seems to converge), but for that you have the integrator configs to change its behaviour:

See the tutorial how to change the config: ROOT: tutorials/roofit/roofit/rf901_numintconfig.C File Reference
And you could use specifically the 1D config here:

  *** RooIntegrator1D ***
  Capabilities: 1-D 
  Configuration: 
    1)        sumRule = Trapezoid(idx = 0)
    2)  extrapolation = Wynn-Epsilon(idx = 1)
    3)       maxSteps = 20
    4)       minSteps = 999
    5)       fixSteps = 0

to force the minimum number of steps to a higher value. ~~999 looks to be a value meaning “deactivated”, but if you set it to 10 or something, the integrator would be forced to run these number of refinements before returning the value.~~
Update: The minimum number of steps only has an effect in cases where the integrator believes that the result might be zero. You can, however, force it to run longer by setting lower values for epsilon, the required accuracy.
A brute-force method is of course to fix the number of steps to some value of your choosing.

Ako_b · September 17, 2025, 11:35am

Thanks for help, @StephanH !

In meanwhile, I’ve constructed pathological example with periodic function to show that aliasing effect can hit RooFit integrators too.

roofit_test_create_projection_aliasing.cxx (2.0 KB)