Home | News | Documentation | Download

Normalization integrals in a 3D conditional pdf

Dear all,

I hit a wall when trying to fit data with a 3D pdf of the form f(x)*g(y)*h(z|(x,y)). It seemed like the fit simply got stuck (burning way too many CPU hours). Attached, you find a simple MWE derived from that fit which I hope reproduces the problem. In the end I think it boils down to: “Can you please have a look at my model and tell me how to set it up correctly” - sorry for that in advance!

In the example, 2D histograms from the 3D model are created and plotted. This works in (x,z), (y,z), but not in (x,y) - or simply takes too long.
I’ll add two comments here to highlight parts of the code I don’t understand or am unsure about.

  • In line 36 I used RooFit::Conditional(*hz,RooArgSet(x,y),true) but don’t understand what true does. I can see it’s not the default argument and it eventually allows to draw the (x,z) and (y,z) projections, which would otherwise be stuck as well. You can also see below that the stdout is different in both cases.
  • The factory command "PROD::factory_model(hx,hy,hz|{x,y})" (which I also used for the failed fit described in the beginning) does not seem to work as expected, as seen by printing the RooWorkspace

For what it’s worth, let me add the printout when setting the flag described in the first point to true/false. This raised suspicion for me that it is related to calculating normalization integrals.

depsAreCond flag off (default), (x,z) projection

[#1] INFO:Caching -- RooAbsCachedPdf::getCache(hx) creating new cache 0x3959520 with pdf hx_2_CONV_hx_3_CACHE_Obs[x] for nset () with code 0 from preexisting content.
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(hz) creating new cache 0x3935b20 with pdf hz_2_CONV_hz_3_CACHE_Obs[z] for nset () with code 0 from preexisting content.
[#1] INFO:NumericIntegration -- RooRealIntegral::init([hx_X_hz]_Norm[x]_denominator_Int[x]) using numeric integrator RooIntegrator1D to calculate Int(x)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(hz_2_Int[z]) using numeric integrator RooIntegrator1D to calculate Int(z)

depsAreCond flag on, (x,z) projection

[#1] INFO:Caching -- RooAbsCachedPdf::getCache(hx) creating new cache 0x4758c50 with pdf hx_2_CONV_hx_3_CACHE_Obs[x] for nset () with code 0 from preexisting content.
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(hx) creating new cache 0x4744cf0 with pdf hx_2_CONV_hx_3_CACHE_Obs[x] for nset (x) with code 1 from preexisting content.
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(hz) creating new cache 0x47a0a50 with pdf hz_2_CONV_hz_3_CACHE_Obs[z] for nset () with code 0 from preexisting content.
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(hz) creating new cache 0x47c58b0 with pdf hz_2_CONV_hz_3_CACHE_Obs[z] for nset (z) with code 1 from preexisting content.
[#1] INFO:NumericIntegration -- RooRealIntegral::init(hz_2_Int[z]) using numeric integrator RooIntegrator1D to calculate Int(z)

and finally (x,y) which always freezes

[#1] INFO:Caching -- RooAbsCachedPdf::getCache(hx) creating new cache 0x47d4740 with pdf hx_2_CONV_hx_3_CACHE_Obs[x] for nset () with code 0 from preexisting content.
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(hz) creating new cache 0x47e0ed0 with pdf hz_2_CONV_hz_3_CACHE_Obs[z] for nset () with code 0 from preexisting content.
[#1] INFO:NumericIntegration -- RooRealIntegral::init(hz_2_Int[y,z]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(y,z)
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(hz) creating new cache 0x48e39b0 with pdf hz_2_CONV_hz_3_CACHE_Obs[y,z] for nset (y) with code 2
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(hy) creating new cache 0x565e440 with pdf hy_2_CONV_hy_3_CACHE_Obs[y] for nset () with code 0 from preexisting content.
[#1] INFO:NumericIntegration -- RooRealIntegral::init([hx_X_hz]_Norm[x,y]_denominator_Int[x]) using numeric integrator RooIntegrator1D to calculate Int(x)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(hz_2_Int[z]) using numeric integrator RooIntegrator1D to calculate Int(z)
[#1] INFO:NumericIntegration -- RooRealIntegral::init([hy_X_hz]_Norm[x,y]_denominator_Int[x,y]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(y,x)

Thanks in advance for any help with this!
Marian

conditional_pdf_toy.cpp (4.2 KB)

May be @moneta or @StephanH can help you.

Hello,

What I can say for now is that it goes into an n-dimensional integrator and doesn’t come back.
That’s the one:
https://root.cern/doc/master/classROOT_1_1Math_1_1AdaptiveIntegratorMultiDim.html

We still have to figure out which function is actually being integrated here. What I see is that the Fourier convolution is part of it, so that’s a heavy load because a numerical integration algorithm is evaluating a numerical convolution algorithm.
Still, eventually, the integration should converge. Is function “ugly” in the sense that it is strongly peaked or something like that? It could be that changing the integration strategy can solve this, but I have to look up which integrators could help here.

If you want to play with integrators, have a look at rf901_numintconfig.C.

Regarding

The boolean inverts the “conditional”:
Conditional(pdf, x, true) = pdf is conditional on x. That is, don’t integrate it over x to normalise.
Conditional(pdf, x, false) = pdf is conditional on all but x (default)

In your case, I haven’t figured out, yet, why this requires integrating hy*hz over xy to plot in xy, but it takes forever because the numerical integrator will have to run the numerical convolutions for each step (and there’s lots of steps to sample the xy space).

Right now, I’m trying the RooMCIntegrator. It’s brute-force and also slow, but maybe it will work.

  RooNumIntConfig* intConfig = RooAbsReal::defaultIntegratorConfig();
  intConfig->method2D().setLabel("RooMCIntegrator");

Hello Stephan,

Thank’s for the hint to configure integrators! I’ll have a look at this later.

Compared to the given ranges, the function is not strongly peaked. I wanted to use the model for a mass fit where the resolutions correspond to about 1/5 of each of the ranges.
It is true that the Fourier convolutions are quite heavy. Nonetheless I was surprised that they (in my case) tend to be more stable than the raw Hypatia. I’m using the convolution to follow the recipe given in sections 5 and 6 of the Hypatia paper.
By the way: it’s great to see (an improved version of) the Hypatia function in RooFit! Thanks a lot for this :slight_smile:

Thanks also for your second post! It is clearer to me now what happens, although I have to sit down and play a bit with it to understand more.

Cheers,
Marian

Hi again,

after thinking about it for a bit, the AdaptiveIntegrator is probably still a good bet, but the default is to run until it reaches 1 000 000 function calls or so. Given that a couple of thousand already take a minute, that’s obviously not feasible. You can change this setting, but you might have to sacrifice a bit of accuracy.
I will talk to @moneta to check if early stopping of that integrator is safe.

EDIT
I have some more on the early stopping. The adaptive integrator is supposed to run for either 100000 steps or until a relative error < 1.E-7 is reached.
After 2000 steps (about 5 mins) I see a precision of ~5.E-5. So this might work:

  RooNumIntConfig* intConfig = RooAbsReal::defaultIntegratorConfig();
  intConfig->setEpsRel(1.E-5);
  intConfig->getConfigSection("RooAdaptiveIntegratorND").setRealValue("maxEval2D", 5000);

With this, you could get acceptable precision and run times.

Hi Stephan,

These settings are very helpful for debugging! The machine I’m using doesn’t seem to be that powerful and the 5 minutes you quote are 15 for me :slight_smile: However, I also got useful results with setting maxEval2D to 100 (~ 2 minutes).
A few questions on the side: can the normalization run parallelized? Which settings do you use for debugging, e.g. to see the number of steps (is it just a matter of increasing the verbosity)? How often is the normalization done in a fit? Can it be cached somehow?

Running the toy with these settings revealed some more practical problems: the model does not look like expected - independent of the settings in RooFit::Conditional or the integrator config. I tried both settings of the depsAreCond flag, and without RooFit::Conditional (i.e. integrate all observables, if i understood correctly). depsAreCond=false did something extremely weird in the (x,z), (y,z) projections and I didn’t even look at (x,y). Setting depsAreCond=true or not using RooFit::Conditional at all gave the same results: the (x,z), (y,z) projections look like expected (in the projections I’m plotting you might call it an ellipse with sloped major axis). For the (x,y) projection, I was expecting something similar to a bivariate Gaussian centered around 2286 and 1970, but the plot shows that the distribution is rotated and centered at higher x and y values (even outside the range of mu_x and mu_y; fixing everything the effective mean in z depends on, i.e. mus and rhos, doesn’t have an effect). Also, I’m failing to see how you can end up with that (x,y) projection by looking at (x,z) and (y,z).
It is quite likely that I just made a stupid mistake somewhere in putting the model together, but I couldn’t find it by staring at it for a bit. I will try to dissect (reduce dimensions, remove complications like the convolution etc.) it later today.

One last thing is the factory command, i.e. that "PROD::factory_model(hx,hy,hz|{x,y})" produces RooProdPdf::factory_model[ hx * hy * y * hz|x ]. Did you have a chance to look why this happens? My idea to write it like this comes from the rf511 tutorial.
model_xy_hist__x_y.pdf (238.4 KB)
model_xz_hist__x_z.pdf (210.2 KB)
model_yz_hist__y_z.pdf (201.9 KB)

Cheers,
Marian

Hi marian,

Good to hear that you could make some progress. By now, we identified a problem with your setup:
For running a 2D integral, RooFit seems to try to integrate out one dimension. So a 2D numerical integration is running a 1D integration for every point in the 2D space. That will take forever. I will try to get in touch with Wouter Verkerke to ask if there’s a way around this, as running the 2D straight would be faster.

  • Parallel? Not that I knew. If I had more time, I would look into that …
  • I debugged by just adding a cout to the integrator. The benefit of having ROOT ready to compile. :slight_smile:
  • The normalisation is done whenever a parameter changes, so usually in every fit step. It’s actually being cached, but when a parameter changes, one needs to recompute. If partial integrals are used somewhere, the caches don’t get wiped if they don’t depend on the parameter that changed. In this respect, RooFit usually does the right thing except for one part that’s also on my list:
    https://sft.its.cern.ch/jira/browse/ROOT-10519
    The bad news: When using an FFT convolution, this actually affects you … :roll_eyes:

Concerning the model, I don’t know how to help at the moment. Maybe you try again to describe in different words what the model is supposed to do, and that might spark an idea.

I don’t see the |{x,y} syntax in the tutorial. It’s probably in rf512, but here you see that the syntax is not intuitive, because the conditional PDF comes after the conditional observable. The result you report therefore is actually what I would expect after looking at the tutorial. Defining the PDF as conditional on two observables looks like a thing that has to be done manually (I doubt that PROD::factory_model(hx,hy |{x,y},hz) can work, but you can actually invert the conditional using PROD::factory_model(hx,hy |~z,hz). That should mean something like conditional on all but z, i.e. x,y).

Hi Stephan,

Thanks for the explanation and pointers!

I played a bit with the factory-code and found something that produced the plots I had in mind (see (x,y) below, the other projections are the same as before). The code that eventually worked is "PROD::factory_model(hx,hy,hz|~z)", which makes sense from what you explained (conditional on all but z).

Now, also the changes to RooNumIntConfig don’t play a role any longer. It seems that only a 1D numeric integral in z is carried out for the (x,z) and (y,z) projections, while (x,y) pops put immediately.

Let’s see how the model does in a fit now :slight_smile:

Cheers,
Marian

model_xy_hist__x_y.pdf (244.4 KB)

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.