Action #1303

Improve fitting convergence behavior for shell models

Added by Knödlseder Jürgen almost 10 years ago. Updated over 9 years ago.

Status:ClosedStart date:07/28/2014
Priority:NormalDue date:
Assigned To:Knödlseder Jürgen% Done:

100%

Category:-
Target version:00-09-00
Duration:

Description

The shell model fit requires a large number of iterations, indicating problem with convergence which are likely related to numerical noise in the computation of the spatial parameter gradients (see #1299). The convergence behavior of the model should be improved.

Recurrence

No recurrence.

Related issues

Related to GammaLib - Bug #1299: Fitting problems with radial disk model Closed 07/25/2014

History

#1 Updated by Knödlseder Jürgen over 9 years ago

  • Status changed from New to In Progress
  • Assigned To set to Knödlseder Jürgen

I studied the impact of the transition points in the integration scheme on the result. It turned out that the transition point related to the shell radius has not been set correctly (the shell radius was retrieved in degrees which the method expected a value in radians).

Originally, the was a single active transition point at transition_point = delta_max - zeta. This gave:

>Iteration   0: -logL=35267.348, Lambda=1.0e-03
  Parameter "Width" drives optimization step (step=0.160656)
  Parameter "Width" hits minimum: 0.01 < 0.01 (1)
 Iteration   1: -logL=35267.348, Lambda=1.0e-03, delta=-46.510, max(|grad|)=16173.840440 [Width:3] (stalled)
  Parameter "Width" does not drive optimization step anymore.
  Parameter "Width" hits minimum: -0.419089 < 0.01 (2)
 Iteration   2: -logL=35267.348, Lambda=1.0e-02, delta=-766.211, max(|grad|)=4695.710231 [Width:3] (stalled)
  Parameter "Width" drives optimization step (step=0.29825)
  Parameter "Width" hits minimum: 0.01 < 0.01 (3)
 Iteration   3: -logL=35267.348, Lambda=1.0e-01, delta=-46.387, max(|grad|)=16942.638241 [Width:3] (stalled)
  Parameter "Width" does not drive optimization step anymore.
 Iteration   4: -logL=35267.348, Lambda=1.0e+00, delta=-29.779, max(|grad|)=-4648.801580 [Width:3] (stalled)
>Iteration   5: -logL=35266.887, Lambda=1.0e+01, delta=0.461, max(|grad|)=-1018.887844 [Width:3]
>Iteration   6: -logL=35260.439, Lambda=1.0e+00, delta=6.448, max(|grad|)=428.821260 [Radius:2]
>Iteration   7: -logL=35258.384, Lambda=1.0e-01, delta=2.055, max(|grad|)=349.658308 [Radius:2]
 Iteration   8: -logL=35258.384, Lambda=1.0e-02, delta=-4.818, max(|grad|)=130.663905 [Radius:2] (stalled)
 Iteration   9: -logL=35258.384, Lambda=1.0e-01, delta=-8.284, max(|grad|)=-445.029473 [Width:3] (stalled)
 Iteration  10: -logL=35258.384, Lambda=1.0e+00, delta=-0.831, max(|grad|)=189.062455 [Radius:2] (stalled)
>Iteration  11: -logL=35258.159, Lambda=1.0e+01, delta=0.225, max(|grad|)=213.675143 [Radius:2]
 Iteration  12: -logL=35258.159, Lambda=1.0e+00, delta=-1.797, max(|grad|)=2028.395494 [Width:3] (stalled)
>Iteration  13: -logL=35258.089, Lambda=1.0e+01, delta=0.069, max(|grad|)=137.837075 [Radius:2]
 Iteration  14: -logL=35258.089, Lambda=1.0e+00, delta=-0.125, max(|grad|)=514.946557 [Width:3] (stalled)
>Iteration  15: -logL=35258.057, Lambda=1.0e+01, delta=0.032, max(|grad|)=86.939405 [Radius:2]
>Iteration  16: -logL=35257.961, Lambda=1.0e+00, delta=0.096, max(|grad|)=78.540429 [Radius:2]
 Iteration  17: -logL=35257.961, Lambda=1.0e-01, delta=-1.318, max(|grad|)=-243.413734 [Width:3] (stalled)
 Iteration  18: -logL=35257.961, Lambda=1.0e+00, delta=-0.075, max(|grad|)=197.304107 [Width:3] (stalled)
>Iteration  19: -logL=35257.950, Lambda=1.0e+01, delta=0.011, max(|grad|)=32.384030 [Radius:2]
 Iteration  20: -logL=35257.950, Lambda=1.0e+00, delta=-0.015, max(|grad|)=125.138318 [Width:3] (stalled)
 Iteration  21: -logL=35257.950, Lambda=1.0e+01, delta=-0.001, max(|grad|)=16.279204 [Radius:2] (stalled)
 Iteration  22: -logL=35257.950, Lambda=1.0e+02, delta=-0.000, max(|grad|)=29.952981 [Radius:2] (stalled)
 Iteration  23: -logL=35257.950, Lambda=1.0e+03, delta=-0.000, max(|grad|)=32.154405 [Radius:2] (stalled)
 Iteration  24: -logL=35257.950, Lambda=1.0e+04, delta=-0.000, max(|grad|)=32.361191 [Radius:2] (stalled)
 Iteration  25: -logL=35257.950, Lambda=1.0e+05, delta=-0.000, max(|grad|)=32.381748 [Radius:2] (stalled)
 Iteration  26: -logL=35257.950, Lambda=1.0e+06, delta=-0.000, max(|grad|)=32.381526 [Radius:2] (stalled)
 Iteration  27: -logL=35257.950, Lambda=1.0e+07, delta=-0.000, max(|grad|)=32.381506 [Radius:2] (stalled)
 Iteration  28: -logL=35257.950, Lambda=1.0e+08, delta=-0.000, max(|grad|)=32.381501 [Radius:2] (stalled)
 Iteration  29: -logL=35257.950, Lambda=1.0e+09, delta=-0.000, max(|grad|)=32.381504 [Radius:2] (stalled)
=== GOptimizerLM ===
 Optimized function value ..: 35257.950
 Absolute precision ........: 0.005
 Optimization status .......: stalled
 Number of parameters ......: 13
 Number of free parameters .: 9
 Number of iterations ......: 29
 Lambda ....................: 1e+10
 Number of models ..........: 2
 Number of parameters ......: 13
=== GModelSky ===
 Name ......................: Gaussian Crab
 Instruments ...............: all
 Instrument scale factors ..: unity
 Observation identifiers ...: all
 Model type ................: ExtendedSource
 Model components ..........: "ShellFunction" * "PowerLaw" * "Constant" 
 Number of parameters ......: 8
 Number of spatial par's ...: 4
  RA .......................: 83.6324 +/- 0.00551318 [-360,360] deg (free,scale=1)
  DEC ......................: 22.0203 +/- 0.004991 [-90,90] deg (free,scale=1)
  Radius ...................: 0.280905 +/- 0.00712947 [0.01,10] deg (free,scale=1)
  Width ....................: 0.125087 +/- 0.00971526 [0.01,10] deg (free,scale=1)
 Number of spectral par's ..: 3
  Prefactor ................: 5.84407e-16 +/- 2.19786e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
  Index ....................: -2.44876 +/- 0.026312 [-0,-5]  (free,scale=-1,gradient)
  PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
 Number of temporal par's ..: 1
  Constant .................: 1 (relative value) (fixed,scale=1,gradient)
=== GCTAModelRadialAcceptance ===
 Name ......................: Background
 Instruments ...............: CTA
 Instrument scale factors ..: unity
 Observation identifiers ...: all
 Model type ................: "Gaussian" * "PowerLaw" * "Constant" 
 Number of parameters ......: 5
 Number of radial par's ....: 1
  Sigma ....................: 3.03552 +/- 0.0398728 [0.01,10] deg2 (free,scale=1,gradient)
 Number of spectral par's ..: 3
  Prefactor ................: 6.00756e-05 +/- 1.79609e-06 [0,0.001] ph/cm2/s/MeV (free,scale=1e-06,gradient)
  Index ....................: -1.85691 +/- 0.0166786 [-0,-5]  (free,scale=-1,gradient)
  PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
 Number of temporal par's ..: 1
  Constant .................: 1 (relative value) (fixed,scale=1,gradient)
 Elapsed time ..............: 115.498 sec

Adding a transition point at transition_point = src_max - delta_max that takes care of the case that the Psf is fully container in the model gave:

>Iteration   0: -logL=35266.697, Lambda=1.0e-03
  Parameter "Width" drives optimization step (step=0.220368)
  Parameter "Width" hits minimum: 0.01 < 0.01 (1)
 Iteration   1: -logL=35266.697, Lambda=1.0e-03, delta=-37.338, max(|grad|)=18086.039507 [Width:3] (stalled)
  Parameter "Width" does not drive optimization step anymore.
  Parameter "Width" hits minimum: -0.283135 < 0.01 (2)
 Iteration   2: -logL=35266.697, Lambda=1.0e-02, delta=-598.127, max(|grad|)=8004.340348 [Width:3] (stalled)
  Parameter "Width" drives optimization step (step=0.377765)
  Parameter "Width" hits minimum: 0.01 < 0.01 (3)
 Iteration   3: -logL=35266.697, Lambda=1.0e-01, delta=-36.883, max(|grad|)=19091.473377 [Width:3] (stalled)
  Parameter "Width" does not drive optimization step anymore.
 Iteration   4: -logL=35266.697, Lambda=1.0e+00, delta=-4.235, max(|grad|)=1228.332795 [Width:3] (stalled)
>Iteration   5: -logL=35265.425, Lambda=1.0e+01, delta=1.272, max(|grad|)=-778.103446 [Width:3]
>Iteration   6: -logL=35260.036, Lambda=1.0e+00, delta=5.389, max(|grad|)=404.556811 [Radius:2]
>Iteration   7: -logL=35259.115, Lambda=1.0e-01, delta=0.922, max(|grad|)=215.402269 [Width:3]
 Iteration   8: -logL=35259.115, Lambda=1.0e-02, delta=-1.206, max(|grad|)=-708.637607 [Width:3] (stalled)
 Iteration   9: -logL=35259.115, Lambda=1.0e-01, delta=-4.082, max(|grad|)=-658.710119 [Width:3] (stalled)
>Iteration  10: -logL=35259.102, Lambda=1.0e+00, delta=0.012, max(|grad|)=-163.582591 [Width:3]
 Iteration  11: -logL=35259.102, Lambda=1.0e-01, delta=-3.654, max(|grad|)=-283.420540 [Width:3] (stalled)
>Iteration  12: -logL=35258.933, Lambda=1.0e+00, delta=0.169, max(|grad|)=110.089966 [Width:3]
 Iteration  13: -logL=35258.933, Lambda=1.0e-01, delta=-1.349, max(|grad|)=274.215424 [Width:3] (stalled)
>Iteration  14: -logL=35258.895, Lambda=1.0e+00, delta=0.038, max(|grad|)=-74.399458 [Width:3]
 Iteration  15: -logL=35258.895, Lambda=1.0e-01, delta=-8.898, max(|grad|)=-662.409464 [Width:3] (stalled)
 Iteration  16: -logL=35258.895, Lambda=1.0e+00, delta=-0.004, max(|grad|)=36.818054 [Width:3] (stalled)
>Iteration  17: -logL=35258.876, Lambda=1.0e+01, delta=0.019, max(|grad|)=-58.062935 [Width:3]
 Iteration  18: -logL=35258.876, Lambda=1.0e+00, delta=-0.008, max(|grad|)=45.512503 [Width:3] (stalled)
>Iteration  19: -logL=35258.868, Lambda=1.0e+01, delta=0.008, max(|grad|)=-46.287111 [Width:3]
 Iteration  20: -logL=35258.868, Lambda=1.0e+00, delta=-0.005, max(|grad|)=47.328912 [Width:3] (stalled)
>Iteration  21: -logL=35258.859, Lambda=1.0e+01, delta=0.009, max(|grad|)=-42.130132 [Width:3]
 Iteration  22: -logL=35258.859, Lambda=1.0e+00, delta=-0.020, max(|grad|)=52.015101 [Width:3] (stalled)
>Iteration  23: -logL=35258.854, Lambda=1.0e+01, delta=0.005, max(|grad|)=-22.348328 [Width:3]
 Iteration  24: -logL=35258.854, Lambda=1.0e+00, delta=-0.004, max(|grad|)=53.393351 [Width:3] (stalled)
>Iteration  25: -logL=35258.852, Lambda=1.0e+01, delta=0.002, max(|grad|)=8.995648 [Radius:2]
=== GOptimizerLM ===
 Optimized function value ..: 35258.852
 Absolute precision ........: 0.005
 Optimization status .......: converged
 Number of parameters ......: 13
 Number of free parameters .: 9
 Number of iterations ......: 25
 Lambda ....................: 1
=== GModels ===
 Number of models ..........: 2
 Number of parameters ......: 13
=== GModelSky ===
 Name ......................: Gaussian Crab
 Instruments ...............: all
 Instrument scale factors ..: unity
 Observation identifiers ...: all
 Model type ................: ExtendedSource
 Model components ..........: "ShellFunction" * "PowerLaw" * "Constant" 
 Number of parameters ......: 8
 Number of spatial par's ...: 4
  RA .......................: 83.632 +/- 0.00570752 [-360,360] deg (free,scale=1)
  DEC ......................: 22.0193 +/- 0.00513503 [-90,90] deg (free,scale=1)
  Radius ...................: 0.274423 +/- 0.0100412 [0.01,10] deg (free,scale=1)
  Width ....................: 0.133876 +/- 0.0136274 [0.01,10] deg (free,scale=1)
 Number of spectral par's ..: 3
  Prefactor ................: 5.86297e-16 +/- 2.20807e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
  Index ....................: -2.44825 +/- 0.0262896 [-0,-5]  (free,scale=-1,gradient)
  PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
 Number of temporal par's ..: 1
  Constant .................: 1 (relative value) (fixed,scale=1,gradient)
=== GCTAModelRadialAcceptance ===
 Name ......................: Background
 Instruments ...............: CTA
 Instrument scale factors ..: unity
 Observation identifiers ...: all
 Model type ................: "Gaussian" * "PowerLaw" * "Constant" 
 Number of parameters ......: 5
 Number of radial par's ....: 1
  Sigma ....................: 3.03743 +/- 0.0399576 [0.01,10] deg2 (free,scale=1,gradient)
 Number of spectral par's ..: 3
  Prefactor ................: 5.9936e-05 +/- 1.79318e-06 [0,0.001] ph/cm2/s/MeV (free,scale=1e-06,gradient)
  Index ....................: -1.85782 +/- 0.0166832 [-0,-5]  (free,scale=-1,gradient)
  PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
 Number of temporal par's ..: 1
  Constant .................: 1 (relative value) (fixed,scale=1,gradient)
 Elapsed time ..............: 122.346 sec

After correcting the transition point related to the shell radius I got the following:

>Iteration   0: -logL=35265.099, Lambda=1.0e-03
>Iteration   1: -logL=35261.796, Lambda=1.0e-03, delta=3.303, max(|grad|)=45.073735 [RA:0]
>Iteration   2: -logL=35261.742, Lambda=1.0e-04, delta=0.054, max(|grad|)=-12.689001 [RA:0]
>Iteration   3: -logL=35261.739, Lambda=1.0e-05, delta=0.003, max(|grad|)=3.982421 [RA:0]
=== GOptimizerLM ===
 Optimized function value ..: 35261.739
 Absolute precision ........: 0.005
 Optimization status .......: converged
 Number of parameters ......: 13
 Number of free parameters .: 9
 Number of iterations ......: 3
 Lambda ....................: 1e-06
=== GModels ===
 Number of models ..........: 2
 Number of parameters ......: 13
=== GModelSky ===
 Name ......................: Gaussian Crab
 Instruments ...............: all
 Instrument scale factors ..: unity
 Observation identifiers ...: all
 Model type ................: ExtendedSource
 Model components ..........: "ShellFunction" * "PowerLaw" * "Constant" 
 Number of parameters ......: 8
 Number of spatial par's ...: 4
  RA .......................: 83.6326 +/- 0.00538575 [-360,360] deg (free,scale=1)
  DEC ......................: 22.0195 +/- 0.00484489 [-90,90] deg (free,scale=1)
  Radius ...................: 0.285584 +/- 0.0125695 [0.01,10] deg (free,scale=1)
  Width ....................: 0.116752 +/- 0.0182694 [0.05,10] deg (free,scale=1)
 Number of spectral par's ..: 3
  Prefactor ................: 5.81654e-16 +/- 2.19191e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
  Index ....................: -2.44849 +/- 0.0263003 [-0,-5]  (free,scale=-1,gradient)
  PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
 Number of temporal par's ..: 1
  Constant .................: 1 (relative value) (fixed,scale=1,gradient)
=== GCTAModelRadialAcceptance ===
 Name ......................: Background
 Instruments ...............: CTA
 Instrument scale factors ..: unity
 Observation identifiers ...: all
 Model type ................: "Gaussian" * "PowerLaw" * "Constant" 
 Number of parameters ......: 5
 Number of radial par's ....: 1
  Sigma ....................: 3.03832 +/- 0.0399827 [0.01,10] deg2 (free,scale=1,gradient)
 Number of spectral par's ..: 3
  Prefactor ................: 5.98494e-05 +/- 1.78956e-06 [0,0.001] ph/cm2/s/MeV (free,scale=1e-06,gradient)
  Index ....................: -1.85884 +/- 0.0166747 [-0,-5]  (free,scale=-1,gradient)
  PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
 Number of temporal par's ..: 1
  Constant .................: 1 (relative value) (fixed,scale=1,gradient)
 Elapsed time ..............: 40.805 sec

This is perfect! So it was just a bug in the transition point computation.

#2 Updated by Knödlseder Jürgen over 9 years ago

  • Status changed from In Progress to Closed
  • Target version set to 00-09-00
  • % Done changed from 0 to 100
  • Remaining (hours) set to 0.0

Also available in: Atom PDF