Bug #1496
Fit errors in log parabola fit are too large
| Status: | Rejected | Start date: | 07/01/2015 | |
|---|---|---|---|---|
| Priority: | Normal | Due date: | ||
| Assigned To: |  Knödlseder Jürgen | % Done: | 0% | |
| Category: | - | |||
| Target version: | 1.0.0 | |||
| Duration: |
Description
When fitting all parameters, the fit errors of the log parabola model are too large.
The problem can be reproduced using this XML file: crab_logparabola.xml.
The ctlike run produces the following output:
2015-07-01T06:03:36: +=================================+ 2015-07-01T06:03:36: | Maximum likelihood optimisation | 2015-07-01T06:03:36: +=================================+ 2015-07-01T06:03:36: >Iteration 0: -logL=133492.928, Lambda=1.0e-03 2015-07-01T06:03:36: >Iteration 1: -logL=133489.599, Lambda=1.0e-03, delta=3.330, max(|grad|)=3.544902 [Index:8] 2015-07-01T06:03:36: >Iteration 2: -logL=133489.598, Lambda=1.0e-04, delta=0.000, max(|grad|)=0.013988 [Index:8] 2015-07-01T06:03:37: 2015-07-01T06:03:37: +=========================================+ 2015-07-01T06:03:37: | Maximum likelihood optimization results | 2015-07-01T06:03:37: +=========================================+ 2015-07-01T06:03:37: === GOptimizerLM === 2015-07-01T06:03:37: Optimized function value ..: 133489.598 2015-07-01T06:03:37: Absolute precision ........: 0.005 2015-07-01T06:03:37: Acceptable value decrease .: 2 2015-07-01T06:03:37: Optimization status .......: converged 2015-07-01T06:03:37: Number of parameters ......: 11 2015-07-01T06:03:37: Number of free parameters .: 6 2015-07-01T06:03:37: Number of iterations ......: 2 2015-07-01T06:03:37: Lambda ....................: 1e-05 2015-07-01T06:03:37: Maximum log likelihood ....: -133489.598 2015-07-01T06:03:37: Observed events (Nobs) ...: 57339.000 2015-07-01T06:03:37: Predicted events (Npred) ..: 57337.998 (Nobs - Npred = 1.00155) 2015-07-01T06:03:37: === GModels === 2015-07-01T06:03:37: Number of models ..........: 2 2015-07-01T06:03:37: Number of parameters ......: 11 2015-07-01T06:03:37: === GModelSky === 2015-07-01T06:03:37: Name ......................: Crab 2015-07-01T06:03:37: Instruments ...............: all 2015-07-01T06:03:37: Instrument scale factors ..: unity 2015-07-01T06:03:37: Observation identifiers ...: all 2015-07-01T06:03:37: Model type ................: PointSource 2015-07-01T06:03:37: Model components ..........: "SkyDirFunction" * "LogParabola" * "Constant" 2015-07-01T06:03:37: Number of parameters ......: 7 2015-07-01T06:03:37: Number of spatial par's ...: 2 2015-07-01T06:03:37: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) 2015-07-01T06:03:37: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) 2015-07-01T06:03:37: Number of spectral par's ..: 4 2015-07-01T06:03:37: Prefactor ................: 5.86089e-16 +/- 4.39833e-11 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) 2015-07-01T06:03:37: Index ....................: -2.32353 +/- 4567.39 [-0,-5] (free,scale=-1,gradient) 2015-07-01T06:03:37: Curvature ................: -0.0707068 +/- 0.00274528 [-0.01,-1000] (free,scale=-1,gradient) 2015-07-01T06:03:37: PivotEnergy ..............: 998705 +/- 3.22562e+10 [10000,1e+09] MeV (free,scale=1e+06,gradient) 2015-07-01T06:03:37: Number of temporal par's ..: 1 2015-07-01T06:03:37: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) 2015-07-01T06:03:37: === GCTAModelIrfBackground === 2015-07-01T06:03:37: Name ......................: Background model 2015-07-01T06:03:37: Instruments ...............: CTA 2015-07-01T06:03:37: Instrument scale factors ..: unity 2015-07-01T06:03:37: Observation identifiers ...: all 2015-07-01T06:03:37: Model type ................: "PowerLaw" * "Constant" 2015-07-01T06:03:37: Number of parameters ......: 4 2015-07-01T06:03:37: Number of spectral par's ..: 3 2015-07-01T06:03:37: Prefactor ................: 0.999265 +/- 0.0215484 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient) 2015-07-01T06:03:37: Index ....................: -0.0117382 +/- 0.0123976 [-5,5] (free,scale=1,gradient) 2015-07-01T06:03:37: PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) 2015-07-01T06:03:37: Number of temporal par's ..: 1 2015-07-01T06:03:37: Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
crab_logparabola.xml 
(1.21 KB)
Recurrence
No recurrence.
Related issues
History
#1
Updated by Knödlseder Jürgen over 9 years ago
- File crab_logparabola.xml added
- Description updated (diff)
#2
Updated by Knödlseder Jürgen over 9 years ago
I just disabled the analytical gradient and replaced it by a numerical gradient. This did not solve the issues, meaning that the problem does not come from a wrong analytical gradient computation:
2015-07-01T07:11:37: === GOptimizerLM === 2015-07-01T07:11:37: Optimized function value ..: 133489.598 2015-07-01T07:11:37: Absolute precision ........: 0.005 2015-07-01T07:11:37: Acceptable value decrease .: 2 2015-07-01T07:11:37: Optimization status .......: errors are inaccurate 2015-07-01T07:11:37: Number of parameters ......: 11 2015-07-01T07:11:37: Number of free parameters .: 6 2015-07-01T07:11:37: Number of iterations ......: 2 2015-07-01T07:11:37: Lambda ....................: 1e-05 2015-07-01T07:11:37: Maximum log likelihood ....: -133489.598 2015-07-01T07:11:37: Observed events (Nobs) ...: 57339.000 2015-07-01T07:11:37: Predicted events (Npred) ..: 57337.998 (Nobs - Npred = 1.00154) 2015-07-01T07:11:37: === GModels === 2015-07-01T07:11:37: Number of models ..........: 2 2015-07-01T07:11:37: Number of parameters ......: 11 2015-07-01T07:11:37: === GModelSky === 2015-07-01T07:11:37: Name ......................: Crab 2015-07-01T07:11:37: Instruments ...............: all 2015-07-01T07:11:37: Instrument scale factors ..: unity 2015-07-01T07:11:37: Observation identifiers ...: all 2015-07-01T07:11:37: Model type ................: PointSource 2015-07-01T07:11:37: Model components ..........: "SkyDirFunction" * "LogParabola" * "Constant" 2015-07-01T07:11:37: Number of parameters ......: 7 2015-07-01T07:11:37: Number of spatial par's ...: 2 2015-07-01T07:11:37: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) 2015-07-01T07:11:37: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) 2015-07-01T07:11:37: Number of spectral par's ..: 4 2015-07-01T07:11:37: Prefactor ................: 5.86087e-16 +/- 1.03055e-11 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) 2015-07-01T07:11:37: Index ....................: -2.32353 +/- 1070.16 [-0,-5] (free,scale=-1,gradient) 2015-07-01T07:11:37: Curvature ................: -0.0707068 +/- 0.00274528 [-0.01,-1000] (free,scale=-1,gradient) 2015-07-01T07:11:37: PivotEnergy ..............: 998706 +/- 7.55782e+09 [10000,1e+09] MeV (free,scale=1e+06) 2015-07-01T07:11:37: Number of temporal par's ..: 1 2015-07-01T07:11:37: Normalization ............: 1 (relative value) (fixed,scale=1,gradient) 2015-07-01T07:11:37: === GCTAModelIrfBackground === 2015-07-01T07:11:37: Name ......................: Background model 2015-07-01T07:11:37: Instruments ...............: CTA 2015-07-01T07:11:37: Instrument scale factors ..: unity 2015-07-01T07:11:37: Observation identifiers ...: all 2015-07-01T07:11:37: Model type ................: "PowerLaw" * "Constant" 2015-07-01T07:11:37: Number of parameters ......: 4 2015-07-01T07:11:37: Number of spectral par's ..: 3 2015-07-01T07:11:37: Prefactor ................: 0.999265 +/- 0.0215484 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient) 2015-07-01T07:11:37: Index ....................: -0.0117382 +/- 0.0123976 [-5,5] (free,scale=1,gradient) 2015-07-01T07:11:37: PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient) 2015-07-01T07:11:37: Number of temporal par's ..: 1 2015-07-01T07:11:37: Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
#3
Updated by Knödlseder Jürgen over 9 years ago
I replaced the error computation in ctlike by the Hessian error computation. Surprisingly, this gives basically the same covariance matrix, yet this covariance matrix is not positive definite:
=== GMatrixSparse === (curvature matrix) Number of rows ............: 11 Number of columns .........: 11 Number of nonzero elements : 36 Number of allocated cells .: 548 Memory block size .........: 512 Sparse matrix fill ........: 0.297521 Pending element ...........: 0 Fill stack size ...........: 0 (none) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 0, 0, 1268.55, 5724.77, -14854.5, 16486.8, 0, 43.0204, -79.7296, 0, 0 0, 0, 5724.77, 87060.3, -93910, 65733.2, 0, 467.303, -917.921, 0, 0 0, 0, -14854.5, -93910, 318438, -189252, 0, -917.99, 1841.59, 0, 0 0, 0, 16486.8, 65733.2, -189252, 215500, 0, 520.441, -957.187, 0, 0 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 0, 0, 43.0204, 467.303, -917.99, 520.441, 0, 13279, -21126.9, 0, 0 0, 0, -79.7296, -917.921, 1841.59, -957.187, 0, -21126.9, 40119.7, 0, 0 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 === GMatrixSparse === (Hessian) Number of rows ............: 11 Number of columns .........: 11 Number of nonzero elements : 36 Number of allocated cells .: 512 Memory block size .........: 512 Sparse matrix fill ........: 0.297521 Pending element ...........: 0 Fill stack size ...........: 0 (none) 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 0, 0, 1268.47, 5725.31, -14851.8, 16487.6, 0, 43.0624, -79.2227, 0, 0 0, 0, 5725.31, 87059.9, -93708.9, 65739, 0, 467.462, -917.827, 0, 0 0, 0, -14851.8, -93708.9, 319855, -189275, 0, -917.479, 1846.5, 0, 0 0, 0, 16487.6, 65739, -189275, 215497, 0, 520.326, -958.55, 0, 0 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 0, 0, 43.0624, 467.462, -917.479, 520.326, 0, 13279, -21127, 0, 0 0, 0, -79.2227, -917.827, 1846.5, -958.55, 0, -21127, 40278.9, 0, 0 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 Non-Positive definite hessian matrix encountered. Load diagonal elements with 1e-10. Fit errors may be inaccurate. Non-Positive definite hessian matrix encountered, even after diagonal loading.
#4
Updated by Knödlseder Jürgen over 9 years ago
I played with the PivotEnergy scale factor, and this had an impact on the result.
<parameter name="Scale" scale="1.0" value="1.0e6" min="0.01" max="1.0e9" free="1"/>
gave
2015-07-01T10:38:41: +==================+ 2015-07-01T10:38:41: | Curvature matrix | 2015-07-01T10:38:41: +==================+ 2015-07-01T10:38:41: === GMatrixSparse === 2015-07-01T10:38:41: Number of rows ............: 11 2015-07-01T10:38:41: Number of columns .........: 11 2015-07-01T10:38:41: Number of nonzero elements : 36 2015-07-01T10:38:41: Number of allocated cells .: 548 2015-07-01T10:38:41: Memory block size .........: 512 2015-07-01T10:38:41: Sparse matrix fill ........: 0.297521 2015-07-01T10:38:41: Pending element ...........: 0 2015-07-01T10:38:41: Fill stack size ...........: 0 (none) 2015-07-01T10:38:41: 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 2015-07-01T10:38:41: 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 2015-07-01T10:38:41: 0, 0, 1270.09, 5730.2, -14866.5, 0.0164925, 0, 43.0466, -79.7782, 0, 0 2015-07-01T10:38:41: 0, 0, 5730.2, 87077.9, -93978.4, 0.0657413, 0, 467.369, -918.043, 0, 0 2015-07-01T10:38:41: 0, 0, -14866.5, -93978.4, 318537, -0.189237, 0, -918.235, 1842.07, 0, 0 2015-07-01T10:38:41: 0, 0, 0.0164925, 0.0657413, -0.189237, 2.15387e-07, 0, 0.000520305, -0.000956936, 0, 0 2015-07-01T10:38:41: 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 2015-07-01T10:38:41: 0, 0, 43.0466, 467.369, -918.235, 0.000520305, 0, 13279, -21126.9, 0, 0 2015-07-01T10:38:41: 0, 0, -79.7782, -918.043, 1842.07, -0.000956936, 0, -21126.9, 40119.7, 0, 0 2015-07-01T10:38:41: 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 2015-07-01T10:38:41: 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 2015-07-01T10:38:41: Non-Positive definite curvature matrix encountered. 2015-07-01T10:38:41: Load diagonal elements with 1e-10. Fit errors may be inaccurate. ... 2015-07-01T10:38:41: === GModelSky === 2015-07-01T10:38:41: Name ......................: Crab 2015-07-01T10:38:41: Instruments ...............: all 2015-07-01T10:38:41: Instrument scale factors ..: unity 2015-07-01T10:38:41: Observation identifiers ...: all 2015-07-01T10:38:41: Model type ................: PointSource 2015-07-01T10:38:41: Model components ..........: "SkyDirFunction" * "LogParabola" * "Constant" 2015-07-01T10:38:41: Number of parameters ......: 7 2015-07-01T10:38:41: Number of spatial par's ...: 2 2015-07-01T10:38:41: RA .......................: 83.6331 [-360,360] deg (fixed,scale=1) 2015-07-01T10:38:41: DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1) 2015-07-01T10:38:41: Number of spectral par's ..: 4 2015-07-01T10:38:41: Prefactor ................: 5.85732e-16 +/- 1.36307e-16 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) 2015-07-01T10:38:41: Index ....................: -2.32357 +/- 0.014771 [-0,-5] (free,scale=-1,gradient) 2015-07-01T10:38:41: Curvature ................: -0.0707068 +/- 0.00274528 [-0.01,-1000] (free,scale=-1,gradient) 2015-07-01T10:38:41: PivotEnergy ..............: 998967 +/- 100000 [0.01,1e+09] MeV (free,scale=1,gradient) 2015-07-01T10:38:41: Number of temporal par's ..: 1 2015-07-01T10:38:41: Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
#5
Updated by Knödlseder Jürgen over 9 years ago
Here some more results:
<parameter name="Scale" scale="10.0" value="1.0e5" min="0.01" max="1.0e9" free="1"/> ... 2015-07-01T10:41:09: Prefactor ................: 5.86064e-16 +/- 5.67734e-11 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) 2015-07-01T10:41:09: Index ....................: -2.32353 +/- 5895.79 [-0,-5] (free,scale=-1,gradient) 2015-07-01T10:41:09: Curvature ................: -0.0707068 +/- 0.00274528 [-0.01,-1000] (free,scale=-1,gradient) 2015-07-01T10:41:09: PivotEnergy ..............: 998723 +/- 4.16386e+10 [0.1,1e+10] MeV (free,scale=10,gradient)
<parameter name="Scale" scale="100.0" value="1.0e4" min="0.01" max="1.0e9" free="1"/> ... 2015-07-01T10:42:48: Prefactor ................: 5.86114e-16 +/- 5.04599e-11 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient) 2015-07-01T10:42:48: Index ....................: -2.32353 +/- 5239.71 [-0,-5] (free,scale=-1,gradient) 2015-07-01T10:42:48: Curvature ................: -0.0707068 +/- 0.00274528 [-0.01,-1000] (free,scale=-1,gradient) 2015-07-01T10:42:48: PivotEnergy ..............: 998687 +/- 3.70038e+10 [1,1e+11] MeV (free,scale=100,gradient)
Looks that playing around with the scales does lead to a non positive definite matrix, resulting in diagonal loading which essentially switches off the pivot energy.
#6
Updated by Knödlseder Jürgen over 9 years ago
After some lunch time discussions we concluded that it is impossible to constrain the pivot energy of the log-parabola model as the model parameters are fully degenerated. This issue is thus considered as a non-issue, and can be rejected.
Note, however, that this should be mentioned in the known issues of the ctools documentation.
#7
Updated by Knödlseder Jürgen over 9 years ago
- Status changed from New to Rejected