Problem description¶

Here a number of log-likelihood profiles for the disk parameters RA (top), DEC (middle), and Radius (bottom). The plots from the left to the right show successive zooms with step size 1e-3 (left), 1e-4 (centre) and 1e-5 (right). The red line connects the computed log-likelihood values (red dots), the blue line indicates the gradient obtained by the optimizer. For gradient computation, a step size of 1e-5 has been used. Obviously, the log-likelihood profile is “nice” at a scale of 1e-3, becomes distorted at a scale of 1e-4 and becomes very noisy at a scale of 1e-5.

Reducing the numerical differentiation step-size did not really change anything. So although the profile looks smooth on a large scale, it is not precise enough for accurate numerical derivatives that let the algorithm converge properly.

Npred computation not the problem¶

I checked that the Npred term is not the problem in the unbinned analysis. Below are Npred profiles for a step size of 1e-5 in all three parameters. Everything is smooth. This implies that it must be the probability density computation for the individual events that creates the structure in the log-likelihood profiles.

Below the log-likelihood profiles with a step size of 1e-5 for the reduced relative integration precision of 1e-2. This seems counter intuitive, but with less precision the profiles are smoother!

Binned also better for 1e-2¶

Here now a comparison of the binned analysis, first with the old precision of 1e-5, and then with the new precision of 1e-2. The fit obviously converges much faster and smoother with the reduced integration precision of 1e-2. The results are pretty much the same.

Integration precision 1e-5

>Iteration   0: -logL=18456.210, Lambda=1.0e-03
>Iteration   1: -logL=18452.575, Lambda=1.0e-03, delta=3.635, max(|grad|)=-28.343815 [RA:0]
>Iteration   2: -logL=18452.529, Lambda=1.0e-04, delta=0.045, max(|grad|)=-16.720985 [DEC:1]
>Iteration   3: -logL=18452.529, Lambda=1.0e-05, delta=0.000, max(|grad|)=6.820850 [DEC:1]
>Iteration   4: -logL=18452.528, Lambda=1.0e-06, delta=0.001, max(|grad|)=-2.600544 [RA:0]
 Iteration   5: -logL=18452.528, Lambda=1.0e-07, delta=-0.000, max(|grad|)=3.205072 [RA:0] (stalled)
 Iteration   6: -logL=18452.528, Lambda=1.0e-06, delta=-0.000, max(|grad|)=3.205211 [RA:0] (stalled)
 Iteration   7: -logL=18452.528, Lambda=1.0e-05, delta=-0.000, max(|grad|)=3.205105 [RA:0] (stalled)
 Iteration   8: -logL=18452.528, Lambda=1.0e-04, delta=-0.000, max(|grad|)=3.204754 [RA:0] (stalled)
 Iteration   9: -logL=18452.528, Lambda=1.0e-03, delta=-0.000, max(|grad|)=3.248355 [RA:0] (stalled)
 Iteration  10: -logL=18452.528, Lambda=1.0e-02, delta=-0.000, max(|grad|)=3.240016 [RA:0] (stalled)
 Iteration  11: -logL=18452.528, Lambda=1.0e-01, delta=-0.000, max(|grad|)=3.080144 [RA:0] (stalled)
>Iteration  12: -logL=18452.528, Lambda=1.0e+00, delta=0.000, max(|grad|)=-1.595976 [Radius:2]
>Iteration  13: -logL=18452.528, Lambda=1.0e-01, delta=0.000, max(|grad|)=2.660137 [DEC:1]
 Iteration  14: -logL=18452.528, Lambda=1.0e-02, delta=-0.000, max(|grad|)=0.518520 [DEC:1] (stalled)
 Iteration  15: -logL=18452.528, Lambda=1.0e-01, delta=-0.000, max(|grad|)=2.844791 [RA:0] (stalled)
>Iteration  16: -logL=18452.528, Lambda=1.0e+00, delta=0.000, max(|grad|)=2.905791 [RA:0]
 Iteration  17: -logL=18452.528, Lambda=1.0e-01, delta=-0.000, max(|grad|)=-2.042002 [RA:0] (stalled)
 Iteration  18: -logL=18452.528, Lambda=1.0e+00, delta=-0.000, max(|grad|)=-1.560911 [RA:0] (stalled)
>Iteration  19: -logL=18452.528, Lambda=1.0e+01, delta=0.000, max(|grad|)=-2.206019 [DEC:1]
=== GObservations ===
 Number of observations ....: 1
 Number of predicted events : 3687
=== GModelSky ===
 Number of spatial par's ...: 3
  RA .......................: 83.6324 +/- 0.00390577 [-360,360] deg (free,scale=1)
  DEC ......................: 22.0143 +/- 0.00329028 [-90,90] deg (free,scale=1)
  Radius ...................: 0.201591 +/- 0.00304322 [0.01,10] deg (free,scale=1)
 Number of spectral par's ..: 3
  Prefactor ................: 5.41603e-16 +/- 1.98302e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
  Index ....................: -2.45653 +/- 0.0267726 [-0,-5]  (free,scale=-1,gradient)
=== GCTAModelRadialAcceptance ===
 Number of radial par's ....: 1
  Sigma ....................: 3.15909 +/- 0.0857976 [0.01,10] deg2 (free,scale=1,gradient)
 Number of spectral par's ..: 3
  Prefactor ................: 5.88695e-05 +/- 1.90709e-06 [0,0.001] ph/cm2/s/MeV (free,scale=1e-06,gradient)
  Index ....................: -1.85371 +/- 0.0175121 [-0,-5]  (free,scale=-1,gradient)

Integration precision 1e-2

>Iteration   0: -logL=18456.397, Lambda=1.0e-03
>Iteration   1: -logL=18452.583, Lambda=1.0e-03, delta=3.814, max(|grad|)=-20.005022 [RA:0]
>Iteration   2: -logL=18452.538, Lambda=1.0e-04, delta=0.046, max(|grad|)=5.343857 [RA:0]
>Iteration   3: -logL=18452.537, Lambda=1.0e-05, delta=0.000, max(|grad|)=-1.070713 [RA:0]
>Iteration   4: -logL=18452.537, Lambda=1.0e-06, delta=0.000, max(|grad|)=0.358104 [RA:0]
>Iteration   5: -logL=18452.537, Lambda=1.0e-07, delta=0.000, max(|grad|)=0.263873 [Radius:2]
=== GObservations ===
 Number of observations ....: 1
 Number of predicted events : 3687
=== GModelSky ===
 Number of spatial par's ...: 3
  RA .......................: 83.6325 +/- 0.00390558 [-360,360] deg (free,scale=1)
  DEC ......................: 22.0143 +/- 0.00329027 [-90,90] deg (free,scale=1)
  Radius ...................: 0.201619 +/- 0.00304471 [0.01,10] deg (free,scale=1)
 Number of spectral par's ..: 3
  Prefactor ................: 5.39378e-16 +/- 1.97492e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
  Index ....................: -2.45609 +/- 0.0267721 [-0,-5]  (free,scale=-1,gradient)
=== GCTAModelRadialAcceptance ===
 Number of radial par's ....: 1
  Sigma ....................: 3.15907 +/- 0.0857955 [0.01,10] deg2 (free,scale=1,gradient)
 Number of spectral par's ..: 3
  Prefactor ................: 5.88702e-05 +/- 1.9071e-06 [0,0.001] ph/cm2/s/MeV (free,scale=1e-06,gradient)
  Index ....................: -1.85371 +/- 0.017512 [-0,-5]  (free,scale=-1,gradient)

Gaussian model has problems at all scales!!!¶

Crosschecking then the unbinned Gaussian model as another case I recognized that also here the log-likelihood profiles are bad, and even for an integration precision of 1e-2 have a kink. Here the iterations of the 1e-5 and 1e-2 precision runs. Fitting results were the same for both runs.

Integration precision 1e-5

>Iteration   0: -logL=35065.905, Lambda=1.0e-03
>Iteration   1: -logL=35062.850, Lambda=1.0e-03, delta=3.055, max(|grad|)=26.700169 [Sigma:2]
>Iteration   2: -logL=35062.829, Lambda=1.0e-04, delta=0.022, max(|grad|)=44.862488 [DEC:1]
 Iteration   3: -logL=35062.829, Lambda=1.0e-05, delta=-0.126, max(|grad|)=-39.775759 [Sigma:2] (stalled)
 Iteration   4: -logL=35062.829, Lambda=1.0e-04, delta=-0.126, max(|grad|)=-39.772432 [Sigma:2] (stalled)
 Iteration   5: -logL=35062.829, Lambda=1.0e-03, delta=-0.126, max(|grad|)=-39.687673 [Sigma:2] (stalled)
 Iteration   6: -logL=35062.829, Lambda=1.0e-02, delta=-0.124, max(|grad|)=-39.369992 [Sigma:2] (stalled)
 Iteration   7: -logL=35062.829, Lambda=1.0e-01, delta=-0.104, max(|grad|)=-71.231283 [Sigma:2] (stalled)
 Iteration   8: -logL=35062.829, Lambda=1.0e+00, delta=-0.039, max(|grad|)=-21.992135 [Sigma:2] (stalled)
 Iteration   9: -logL=35062.829, Lambda=1.0e+01, delta=-0.002, max(|grad|)=-6.135290 [Sigma:2] (stalled)
 Iteration  10: -logL=35062.829, Lambda=1.0e+02, delta=-0.000, max(|grad|)=2.564888 [RA:0] (stalled)
 Iteration  11: -logL=35062.829, Lambda=1.0e+03, delta=-0.000, max(|grad|)=43.240077 [Sigma:2] (stalled)
 Iteration  12: -logL=35062.829, Lambda=1.0e+04, delta=-0.000, max(|grad|)=44.858612 [DEC:1] (stalled)

Integration precision 1e-2

>Iteration   0: -logL=35062.491, Lambda=1.0e-03
>Iteration   1: -logL=35059.403, Lambda=1.0e-03, delta=3.087, max(|grad|)=27.067091 [Sigma:2]
>Iteration   2: -logL=35059.383, Lambda=1.0e-04, delta=0.020, max(|grad|)=-2.179918 [Sigma:2]
>Iteration   3: -logL=35059.383, Lambda=1.0e-05, delta=0.000, max(|grad|)=-0.884052 [RA:0]
 Iteration   4: -logL=35059.383, Lambda=1.0e-06, delta=-0.000, max(|grad|)=0.904869 [RA:0] (stalled)
 Iteration   5: -logL=35059.383, Lambda=1.0e-05, delta=-0.000, max(|grad|)=0.904726 [RA:0] (stalled)
 Iteration   6: -logL=35059.383, Lambda=1.0e-04, delta=-0.000, max(|grad|)=0.904932 [RA:0] (stalled)
 Iteration   7: -logL=35059.383, Lambda=1.0e-03, delta=-0.000, max(|grad|)=0.904089 [RA:0] (stalled)
 Iteration   8: -logL=35059.383, Lambda=1.0e-02, delta=-0.000, max(|grad|)=0.893860 [RA:0] (stalled)
 Iteration   9: -logL=35059.383, Lambda=1.0e-01, delta=-0.000, max(|grad|)=0.803572 [RA:0] (stalled)
>Iteration  10: -logL=35059.383, Lambda=1.0e+00, delta=0.000, max(|grad|)=0.253246 [RA:0]
>Iteration  11: -logL=35059.383, Lambda=1.0e-01, delta=0.000, max(|grad|)=-0.104525 [RA:0]
>Iteration  12: -logL=35059.383, Lambda=1.0e-02, delta=0.000, max(|grad|)=0.079574 [RA:0]

And here the log-likelihood profiles, on the left for a precision of 1e-5, on the right for a precision of 1e-2.

Just a note: for the disk model, the azimuthal integration (using the Psf-based system) can be replaced by an analytical model, at least if Aeff is assumed constant over the Psf. Doing this leads to the same behaviour (see also the log-likelihood plot below). This means that the delta integration should be the culprit ...

The following plots indicate why the two integration schemes differ. What is shown below is the evaluation of the response along the radial integration direction for the first event in the test dataset for various integration precisions. For model-based integration this is radially outwards from the model centre (coordinate rho), for the Psf-based integration this is radially outwards from the Psf centre (coordinate delta). Obviously, the shape is not the same, which is likely due to the fact that the Psf integration region is larger than the model region. The kink in the Psf-based integration will obviously lead to problems. One can see how the algorithm refines the integration steps with increasing precision, while for the model-based integration, all precisions are reached with the same step size.

Splitting the integral into two parts¶

I now implemented a split of the integral into two parts, one to the left and another to the right of the kink (which is fully predictable). Using the model system based integration, an integration precision of 1e-5 and a gradient step size of h=2e-4, I obtained the following log-likelihood profiles (for a step size of 1e-5). Left is before split into two part, right is afterwards. Including the split results in a smoother, though not perfect profile.

It is somehow disappointing that the precision problem of the log-likelihood profile cannot be solved. Globally, the Psf system based integration seems even to be worse that the model system based integration. Somehow the model system based integration kernel seems to be smoother, and thus easier to integrate.

Issue solved?¶

Thinking again over it, I’m wondering whether the “noise” may not come from the fact that the integration precision is not guaranteed. All numerical integration methods do some iterative value refinement, and stop when the difference between subsequent iterations falls below a limit. This may lead to different numbers of iterations depending on details of the setup. When now a parameter is varied, the number of iterations may change, leading to a different value at the end.

I checked this idea using the gauss model, unbinned analysis and model system based integration. For a precision of 1e-2, all rho integrations needed 5 iterations, while for a precision of 1e-5 the number of iterations varied between 5 and 7. The calculations were done using the (2,5) noise-robust derivatives for a step-size h=1e-3.

Here the corresponding log-likelihood profiles (from 5 iterations on the left to 7 iterations on the right):

Now the same using the simple difference derivatives for a step-size of h=1e-5.

This looks pretty good. Good a quick convergence even for h=1e-5. Now the same (simple difference derivatives, step-size of h=1e-5) for the disk model:

For shell models there is still an issue to be solved ...

Conclusion: iterative numerical integration methods with uncontrolled iterations depth are not suited for subsequent numerical differentiation. If numerical derivatives are needed, keep the iterations identical.

Still some trouble with Psf system based integration¶

I now checked the same philosophy for the Psf system based integration. I started with the disk model.

It turns out that with Psf system based integration there is still a problem: the log-likelihood profile is now smooth, but there appears a clear step near the minimum in the profile (see left plot below for 5 iterations). This step changes (see right plot for 8 iterations) with an increasing number of iterations for the integration. This could indicate that the step is related to the limited precision of the integration.

As with the change of the parameters the switch point between full and partial containment of the model within the Psf is changing, the function evaluation in the left and right part of the switch point changes, being a source for discontinuities (as every side will have it’s own intrinsic precision).

Testing the Psf system based on a model that has a radius larger than the Psf radius, and consequently that has no switch points, shows that the problem indeed is related to the switch-point: without switch-point the log-likelihood profiles are smooth.

The question then arose whether the Model system based integration has the same problem with the switch point. My conclusion is no, since the switch point for a Model system based integration will always be at the tail of the Psf, i.e. in a region where the integral is small. The switch point occurs when the Psf is not fully container in the Model, and this happens first for the tail. In a Psf system based integration, the situation is different as in a disk model the model is high up to the boundary (and the Phi angle is largest at the boundary). Consequently, the switch point is always close to the maximum.

Conclusion: a Psf system based integration has poorer properties due to the occurrence of the switch point at function maximum.

I now switch-off the Psf system based integration for unbinned and binned analysis.

The shell model still remains a tricky model to integrate, and so far the code is not yet fully satisfactory. Some more investigations are needed.

For the moment let’s stay with the code as is. Here now a check of the Irf computation for all model types and analysis. I show here the computing time needed to estimate the total number of counts in a model, and in parentheses the relative precision achieved. Results are shown for unbinned, binned and stacked analysis. For unbinned, the Npred methods are tested, for binned, the Irf methods are tested. Hence this provides a complete measure of the computing speed (in seconds) and of the precision ((estimated-expected)/expected) of the CTA response computation. The test is done using the script benchmark_irf_computation.py.

Close now since follow-up actions have been identified.