Bug #1561: Elliptical Gauss model shows fit convergence problems - GammaLib - CTA IRAP Project Gateway

Bug #1561

Elliptical Gauss model shows fit convergence problems

Added by Knödlseder Jürgen over 8 years ago. Updated over 6 years ago.

Status:

In Progress

Start date:

10/30/2015

Priority:

Normal

Due date:

Assigned To:

Knödlseder Jürgen

% Done:

50%

Category:

Target version:

Duration:

Description

I found convergence problems with the elliptical Gauss model that explain why the fit is so slow. The following model leads to the problems:

<?xml version="1.0" standalone="no"?>
<source_library title="source library">
  <source name="Crab" type="ExtendedSource">
    <spectrum type="PowerLaw">
      <parameter name="Prefactor" scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
      <parameter name="Index"     scale="-1"    value="2.48" min="0.0"   max="+5.0"   free="1"/>
      <parameter name="Scale"     scale="1e6"   value="0.3"  min="0.01"  max="1000.0" free="0"/>
    </spectrum>
    <spatialModel type="EllipticalGauss">
      <parameter name="RA"          scale="1.0" value="83.6331" min="-360"  max="360" free="1"/>
      <parameter name="DEC"         scale="1.0" value="22.0145" min="-90"   max="90"  free="1"/>
      <parameter name="PA"          scale="1.0" value="45.0"    min="-360"  max="360" free="1"/>
      <parameter name="MinorRadius" scale="1.0" value="0.2"     min="0.001" max="10"  free="1"/>
      <parameter name="MajorRadius" scale="1.0" value="0.4"     min="0.001" max="10"  free="1"/>
    </spatialModel>
  </source>
  <source name="Background" type="CTAIrfBackground" instrument="CTA">
    <spectrum type="PowerLaw">
      <parameter name="Prefactor" scale="1.0" value="1.0" min="1e-3" max="1e3"    free="1"/>
      <parameter name="Index"     scale="1.0" value="0.0" min="-5.0" max="+5.0"   free="1"/>
      <parameter name="Scale"     scale="1e6" value="1.0" min="0.01" max="1000.0" free="0"/>
    </spectrum>
  </source>
</source_library>

The problems can be seen in the result of the pull distribution fit:

2015-08-21T14:00:34: +=================================+
2015-08-21T14:00:34: | Maximum likelihood optimisation |
2015-08-21T14:00:34: +=================================+
2015-08-21T14:02:12:  >Iteration   0: -logL=108962.288, Lambda=1.0e-03
2015-08-21T14:03:40:   Iteration   1: -logL=108962.288, Lambda=1.0e-03, delta=-62.781, max(|grad|)=-1422.695134 [RA:0] (stalled)
2015-08-21T14:05:14:   Iteration   2: -logL=108962.288, Lambda=1.0e-02, delta=-61.650, max(|grad|)=-1409.424285 [RA:0] (stalled)
2015-08-21T14:06:54:   Iteration   3: -logL=108962.288, Lambda=1.0e-01, delta=-52.524, max(|grad|)=-1309.322322 [RA:0] (stalled)
2015-08-21T14:08:31:   Iteration   4: -logL=108962.288, Lambda=1.0e+00, delta=-19.640, max(|grad|)=-833.365056 [RA:0] (stalled)
2015-08-21T14:10:11:   Iteration   5: -logL=108963.142, Lambda=1.0e+01, delta=-0.854, max(|grad|)=-901.230893 [RA:0] (stalled)
2015-08-21T14:11:56:  >Iteration   6: -logL=108962.669, Lambda=1.0e+02, delta=0.473, max(|grad|)=-471.078708 [RA:0]
2015-08-21T14:13:41:   Iteration   7: -logL=108962.790, Lambda=1.0e+01, delta=-0.120, max(|grad|)=1241.129876 [RA:0] (stalled)
2015-08-21T14:15:15:  >Iteration   8: -logL=108961.539, Lambda=1.0e+02, delta=1.250, max(|grad|)=1090.299204 [RA:0]
2015-08-21T14:17:00:  >Iteration   9: -logL=108961.068, Lambda=1.0e+01, delta=0.471, max(|grad|)=-540.297475 [RA:0]
2015-08-21T14:18:36:   Iteration  10: -logL=108961.883, Lambda=1.0e+00, delta=-0.815, max(|grad|)=202.976034 [DEC:1] (stalled)
2015-08-21T14:20:15:  >Iteration  11: -logL=108961.391, Lambda=1.0e+01, delta=0.492, max(|grad|)=188.715977 [RA:0]
2015-08-21T14:21:58:  >Iteration  12: -logL=108958.816, Lambda=1.0e+00, delta=2.575, max(|grad|)=426.349305 [DEC:1]
2015-08-21T14:23:41:   Iteration  13: -logL=108958.816, Lambda=1.0e-01, delta=-0.711, max(|grad|)=224.674851 [RA:0] (stalled)
2015-08-21T14:25:24:  >Iteration  14: -logL=108957.509, Lambda=1.0e+00, delta=1.307, max(|grad|)=-617.376738 [DEC:1]
2015-08-21T14:27:05:   Iteration  15: -logL=108957.509, Lambda=1.0e-01, delta=-11.635, max(|grad|)=565.058180 [DEC:1] (stalled)
2015-08-21T14:28:45:   Iteration  16: -logL=108957.509, Lambda=1.0e+00, delta=-3.688, max(|grad|)=205.103811 [DEC:1] (stalled)
2015-08-21T14:30:32:  >Iteration  17: -logL=108956.518, Lambda=1.0e+01, delta=0.991, max(|grad|)=307.205446 [DEC:1]
2015-08-21T14:32:08:   Iteration  18: -logL=108956.518, Lambda=1.0e+00, delta=-0.556, max(|grad|)=803.938079 [DEC:1] (stalled)
2015-08-21T14:33:56:  >Iteration  19: -logL=108956.507, Lambda=1.0e+01, delta=0.010, max(|grad|)=-289.729848 [DEC:1]
2015-08-21T14:35:31:   Iteration  20: -logL=108956.507, Lambda=1.0e+00, delta=-3.233, max(|grad|)=129.124075 [DEC:1] (stalled)
2015-08-21T14:37:14:   Iteration  21: -logL=108956.507, Lambda=1.0e+01, delta=-0.013, max(|grad|)=380.779002 [DEC:1] (stalled)
2015-08-21T14:38:57:  >Iteration  22: -logL=108956.447, Lambda=1.0e+02, delta=0.060, max(|grad|)=-220.031413 [DEC:1]
2015-08-21T14:40:38:   Iteration  23: -logL=108956.447, Lambda=1.0e+01, delta=-0.005, max(|grad|)=319.798781 [DEC:1] (stalled)
2015-08-21T14:42:17:  >Iteration  24: -logL=108956.409, Lambda=1.0e+02, delta=0.038, max(|grad|)=-162.676580 [DEC:1]
2015-08-21T14:43:56:  >Iteration  25: -logL=108956.403, Lambda=1.0e+01, delta=0.006, max(|grad|)=261.122944 [DEC:1]
2015-08-21T14:45:31:   Iteration  26: -logL=108956.403, Lambda=1.0e+00, delta=-1.164, max(|grad|)=795.467853 [DEC:1] (stalled)
2015-08-21T14:47:05:   Iteration  27: -logL=108956.403, Lambda=1.0e+01, delta=-0.032, max(|grad|)=-302.539921 [DEC:1] (stalled)
2015-08-21T14:48:49:  >Iteration  28: -logL=108956.368, Lambda=1.0e+02, delta=0.035, max(|grad|)=202.644334 [DEC:1]
2015-08-21T14:50:27:   Iteration  29: -logL=108956.368, Lambda=1.0e+01, delta=-0.008, max(|grad|)=-229.737232 [DEC:1] (stalled)
2015-08-21T14:52:12:  >Iteration  30: -logL=108956.347, Lambda=1.0e+02, delta=0.021, max(|grad|)=157.438869 [DEC:1]
2015-08-21T14:53:56:  >Iteration  31: -logL=108956.340, Lambda=1.0e+01, delta=0.007, max(|grad|)=-169.521882 [DEC:1]
2015-08-21T14:55:35:   Iteration  32: -logL=108956.340, Lambda=1.0e+00, delta=-2.594, max(|grad|)=444.524648 [DEC:1] (stalled)
2015-08-21T14:57:10:   Iteration  33: -logL=108956.340, Lambda=1.0e+01, delta=-0.029, max(|grad|)=306.672572 [DEC:1] (stalled)
2015-08-21T14:58:49:  >Iteration  34: -logL=108956.316, Lambda=1.0e+02, delta=0.024, max(|grad|)=-117.537404 [DEC:1]
2015-08-21T15:00:29:   Iteration  35: -logL=108956.316, Lambda=1.0e+01, delta=-0.005, max(|grad|)=236.848867 [DEC:1] (stalled)
2015-08-21T15:02:12:  >Iteration  36: -logL=108956.302, Lambda=1.0e+02, delta=0.014, max(|grad|)=96.053356 [RA:0]
2015-08-21T15:03:53:  >Iteration  37: -logL=108956.293, Lambda=1.0e+01, delta=0.009, max(|grad|)=181.094062 [DEC:1]
2015-08-21T15:05:40:   Iteration  38: -logL=108956.293, Lambda=1.0e+00, delta=-2.074, max(|grad|)=147.701889 [DEC:1] (stalled)
2015-08-21T15:07:15:   Iteration  39: -logL=108956.293, Lambda=1.0e+01, delta=-0.029, max(|grad|)=-251.997590 [DEC:1] (stalled)
2015-08-21T15:08:54:  >Iteration  40: -logL=108956.277, Lambda=1.0e+02, delta=0.016, max(|grad|)=135.134822 [DEC:1]
2015-08-21T15:10:36:   Iteration  41: -logL=108956.277, Lambda=1.0e+01, delta=-0.004, max(|grad|)=-177.703403 [DEC:1] (stalled)
2015-08-21T15:12:14:  >Iteration  42: -logL=108956.267, Lambda=1.0e+02, delta=0.010, max(|grad|)=102.118817 [DEC:1]
2015-08-21T15:13:54:  >Iteration  43: -logL=108956.259, Lambda=1.0e+01, delta=0.009, max(|grad|)=-122.391464 [DEC:1]
2015-08-21T15:15:33:   Iteration  44: -logL=108956.259, Lambda=1.0e+00, delta=-2.033, max(|grad|)=659.108683 [DEC:1] (stalled)
2015-08-21T15:17:17:   Iteration  45: -logL=108956.259, Lambda=1.0e+01, delta=-0.029, max(|grad|)=275.529340 [DEC:1] (stalled)
2015-08-21T15:18:56:  >Iteration  46: -logL=108956.245, Lambda=1.0e+02, delta=0.013, max(|grad|)=87.118714 [RA:0]
2015-08-21T15:20:36:   Iteration  47: -logL=108956.245, Lambda=1.0e+01, delta=-0.004, max(|grad|)=200.456491 [DEC:1] (stalled)
2015-08-21T15:22:17:  >Iteration  48: -logL=108956.238, Lambda=1.0e+02, delta=0.008, max(|grad|)=75.336813 [RA:0]
2015-08-21T15:23:53:  >Iteration  49: -logL=108956.229, Lambda=1.0e+01, delta=0.008, max(|grad|)=145.054215 [DEC:1]
2015-08-21T15:25:37:   Iteration  50: -logL=108956.229, Lambda=1.0e+00, delta=-2.077, max(|grad|)=-254.272515 [DEC:1] (stalled)
2015-08-21T15:27:25:   Iteration  51: -logL=108956.229, Lambda=1.0e+01, delta=-0.029, max(|grad|)=-237.827189 [DEC:1] (stalled)
2015-08-21T15:29:05:  >Iteration  52: -logL=108956.219, Lambda=1.0e+02, delta=0.010, max(|grad|)=104.127190 [DEC:1]
2015-08-21T15:30:43:   Iteration  53: -logL=108956.219, Lambda=1.0e+01, delta=-0.004, max(|grad|)=-158.281731 [DEC:1] (stalled)
2015-08-21T15:32:23:  >Iteration  54: -logL=108956.213, Lambda=1.0e+02, delta=0.006, max(|grad|)=76.317731 [DEC:1]
2015-08-21T15:34:03:  >Iteration  55: -logL=108956.206, Lambda=1.0e+01, delta=0.007, max(|grad|)=-102.711962 [DEC:1]
2015-08-21T15:35:47:   Iteration  56: -logL=108956.206, Lambda=1.0e+00, delta=-1.761, max(|grad|)=753.001479 [DEC:1] (stalled)
2015-08-21T15:37:19:   Iteration  57: -logL=108956.206, Lambda=1.0e+01, delta=-0.032, max(|grad|)=269.654185 [DEC:1] (stalled)
2015-08-21T15:38:59:  >Iteration  58: -logL=108956.196, Lambda=1.0e+02, delta=0.010, max(|grad|)=73.088665 [RA:0]
2015-08-21T15:40:34:   Iteration  59: -logL=108956.196, Lambda=1.0e+01, delta=-0.005, max(|grad|)=185.786016 [DEC:1] (stalled)
2015-08-21T15:42:15:  >Iteration  60: -logL=108956.191, Lambda=1.0e+02, delta=0.005, max(|grad|)=62.420519 [RA:0]
2015-08-21T15:43:57:  >Iteration  61: -logL=108956.185, Lambda=1.0e+01, delta=0.006, max(|grad|)=127.568193 [DEC:1]
2015-08-21T15:45:43:   Iteration  62: -logL=108956.185, Lambda=1.0e+00, delta=-2.007, max(|grad|)=-436.839358 [DEC:1] (stalled)
2015-08-21T15:47:33:   Iteration  63: -logL=108956.185, Lambda=1.0e+01, delta=-0.034, max(|grad|)=-245.600053 [DEC:1] (stalled)
2015-08-21T15:49:12:  >Iteration  64: -logL=108956.176, Lambda=1.0e+02, delta=0.008, max(|grad|)=87.470588 [DEC:1]
2015-08-21T15:50:47:   Iteration  65: -logL=108956.176, Lambda=1.0e+01, delta=-0.006, max(|grad|)=-154.658341 [DEC:1] (stalled)
2015-08-21T15:52:20:  >Iteration  66: -logL=108956.172, Lambda=1.0e+02, delta=0.004, max(|grad|)=61.711204 [DEC:1]
2015-08-21T15:52:20:  
2015-08-21T15:52:20:  +==================+
2015-08-21T15:52:20:  | Curvature matrix |
2015-08-21T15:52:20:  +==================+
2015-08-21T15:52:20:  === GMatrixSparse ===
2015-08-21T15:52:20:   Number of rows ............: 13
2015-08-21T15:52:20:   Number of columns .........: 13
2015-08-21T15:52:20:   Number of nonzero elements : 81
2015-08-21T15:52:20:   Number of allocated cells .: 593
2015-08-21T15:52:20:   Memory block size .........: 512
2015-08-21T15:52:20:   Sparse matrix fill ........: 0.47929
2015-08-21T15:52:20:   Pending element ...........: 0
2015-08-21T15:52:20:   Fill stack size ...........: 0 (none)
2015-08-21T15:52:20:   15839.7, -9820.18, -1.24854, -260.507, 837.681, -0.910593, 19.4204, 0, 0, -10.6466, 26.5364, 0, 0
2015-08-21T15:52:20:   -9820.18, 17709.6, 2.26539, 129.736, 31.9295, -20.1912, -72.4292, 0, 0, -31.9753, 92.9646, 0, 0
2015-08-21T15:52:20:   -1.24854, 2.26539, 0.494084, -0.161682, -1.55698, 0.013647, -0.0163916, 0, 0, -0.065022, 0.0423105, 0, 0
2015-08-21T15:52:20:   -260.507, 129.736, -0.161682, 10731, 1352.88, -244.975, 121.064, 0, 0, 831.063, -879.816, 0, 0
2015-08-21T15:52:20:   837.681, 31.9295, -1.55698, 1352.88, 37436.6, -376.737, 6.83239, 0, 0, 1906.19, -2193.39, 0, 0
2015-08-21T15:52:20:   -0.910593, -20.1912, 0.013647, -244.975, -376.737, 56.3002, -118.254, 0, 0, 175.477, -270.351, 0, 0
2015-08-21T15:52:20:   19.4204, -72.4292, -0.0163916, 121.064, 6.83239, -118.254, 2784.35, 0, 0, 323.287, -1026.12, 0, 0
2015-08-21T15:52:20:   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
2015-08-21T15:52:20:   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
2015-08-21T15:52:20:   -10.6466, -31.9753, -0.065022, 831.063, 1906.19, 175.477, 323.287, 0, 0, 12175.4, -19490.1, 0, 0
2015-08-21T15:52:20:   26.5364, 92.9646, 0.0423105, -879.816, -2193.39, -270.351, -1026.12, 0, 0, -19490.1, 37417.1, 0, 0
2015-08-21T15:52:20:   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
2015-08-21T15:52:20:   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
2015-08-21T15:54:00: 
2015-08-21T15:54:00: +=========================================+
2015-08-21T15:54:00: | Maximum likelihood optimization results |
2015-08-21T15:54:00: +=========================================+
2015-08-21T15:54:00: === GOptimizerLM ===
2015-08-21T15:54:00:  Optimized function value ..: 108956.172
2015-08-21T15:54:00:  Absolute precision ........: 0.005
2015-08-21T15:54:00:  Acceptable value decrease .: 2
2015-08-21T15:54:00:  Optimization status .......: converged
2015-08-21T15:54:00:  Number of parameters ......: 13
2015-08-21T15:54:00:  Number of free parameters .: 9
2015-08-21T15:54:00:  Number of iterations ......: 66
2015-08-21T15:54:00:  Lambda ....................: 10
2015-08-21T15:54:00:  Maximum log likelihood ....: -108956.172
2015-08-21T15:54:00:  Observed events  (Nobs) ...: 16083.000
2015-08-21T15:54:00:  Predicted events (Npred) ..: 16069.565 (Nobs - Npred = 13.4354)
2015-08-21T15:54:00: === GModels ===
2015-08-21T15:54:00:  Number of models ..........: 2
2015-08-21T15:54:00:  Number of parameters ......: 13
2015-08-21T15:54:00: === GModelSky ===
2015-08-21T15:54:00:  Name ......................: Crab
2015-08-21T15:54:00:  Instruments ...............: all
2015-08-21T15:54:00:  Instrument scale factors ..: unity
2015-08-21T15:54:00:  Observation identifiers ...: all
2015-08-21T15:54:00:  Model type ................: ExtendedSource
2015-08-21T15:54:00:  Model components ..........: "EllipticalGauss" * "PowerLaw" * "Constant" 
2015-08-21T15:54:00:  Number of parameters ......: 9
2015-08-21T15:54:00:  Number of spatial par's ...: 5
2015-08-21T15:54:00:   RA .......................: 83.6417 +/- 0.00982085 [-360,360] deg (free,scale=1)
2015-08-21T15:54:00:   DEC ......................: 22.0147 +/- 0.00928538 [-90,90] deg (free,scale=1)
2015-08-21T15:54:00:   PA .......................: 46.0707 +/- 1.42317 [-360,360] deg (free,scale=1)
2015-08-21T15:54:00:   MajorRadius ..............: 0.400028 +/- 0.010374 [0.001,10] deg (free,scale=1)
2015-08-21T15:54:00:   MinorRadius ..............: 0.206587 +/- 0.00546521 [0.001,10] deg (free,scale=1)
2015-08-21T15:54:00:  Number of spectral par's ..: 3
2015-08-21T15:54:00:   Prefactor ................: 5.6221e-16 +/- 1.6123e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
2015-08-21T15:54:00:   Index ....................: -2.51816 +/- 0.0203739 [-0,-5]  (free,scale=-1,gradient)
2015-08-21T15:54:00:   PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
2015-08-21T15:54:00:  Number of temporal par's ..: 1
2015-08-21T15:54:00:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
2015-08-21T15:54:00: === GCTAModelIrfBackground ===
2015-08-21T15:54:00:  Name ......................: Background model
2015-08-21T15:54:00:  Instruments ...............: CTA
2015-08-21T15:54:00:  Instrument scale factors ..: unity
2015-08-21T15:54:00:  Observation identifiers ...: all
2015-08-21T15:54:00:  Model type ................: "PowerLaw" * "Constant" 
2015-08-21T15:54:00:  Number of parameters ......: 4
2015-08-21T15:54:00:  Number of spectral par's ..: 3
2015-08-21T15:54:00:   Prefactor ................: 1.00587 +/- 0.0227185 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient)
2015-08-21T15:54:00:   Index ....................: 0.00228971 +/- 0.012819 [-5,5]  (free,scale=1,gradient)
2015-08-21T15:54:00:   PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
2015-08-21T15:54:00:  Number of temporal par's ..: 1
2015-08-21T15:54:00:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)

Recurrence

No recurrence.

History

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

Status changed from New to In Progress

I could reproduce the problem with a simple ctobssim run followed by ctlike.

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

% Done changed from 0 to 10

I verified that the elliptical Gaussian model is correctly normalized. Note that the normalization only is accurate in the small angle approximation.

I did a special run with a model where the major and minor axis are identical, and hence the position angle is fit. This went well:

2015-10-30T08:40:04: +=================================+
2015-10-30T08:40:04: | Maximum likelihood optimisation |
2015-10-30T08:40:04: +=================================+
2015-10-30T08:40:22:  >Iteration   0: -logL=107816.489, Lambda=1.0e-03
2015-10-30T08:40:45:  >Iteration   1: -logL=107811.316, Lambda=1.0e-03, delta=5.173, max(|grad|)=-39.237446 [MajorRadius:3]
2015-10-30T08:41:08:  >Iteration   2: -logL=107811.263, Lambda=1.0e-04, delta=0.053, max(|grad|)=3.248685 [RA:0]
2015-10-30T08:41:31:  >Iteration   3: -logL=107811.263, Lambda=1.0e-05, delta=0.000, max(|grad|)=-0.348673 [MajorRadius:3]
...
2015-10-30T08:41:54: +=========================================+
2015-10-30T08:41:54: | Maximum likelihood optimization results |
2015-10-30T08:41:54: +=========================================+
2015-10-30T08:41:54: === GOptimizerLM ===
2015-10-30T08:41:54:  Optimized function value ..: 107811.263
2015-10-30T08:41:54:  Absolute precision ........: 0.005
2015-10-30T08:41:54:  Acceptable value decrease .: 2
2015-10-30T08:41:54:  Optimization status .......: converged
2015-10-30T08:41:54:  Number of parameters ......: 13
2015-10-30T08:41:54:  Number of free parameters .: 8
2015-10-30T08:41:54:  Number of iterations ......: 3
2015-10-30T08:41:54:  Lambda ....................: 1e-06
2015-10-30T08:41:54:  Maximum log likelihood ....: -107811.263
2015-10-30T08:41:54:  Observed events  (Nobs) ...: 16051.000
2015-10-30T08:41:54:  Predicted events (Npred) ..: 16050.999 (Nobs - Npred = 0.000532591)
2015-10-30T08:41:54: === GModels ===
2015-10-30T08:41:54:  Number of models ..........: 2
2015-10-30T08:41:54:  Number of parameters ......: 13
2015-10-30T08:41:54: === GModelSky ===
2015-10-30T08:41:54:  Name ......................: Crab
2015-10-30T08:41:54:  Instruments ...............: all
2015-10-30T08:41:54:  Instrument scale factors ..: unity
2015-10-30T08:41:54:  Observation identifiers ...: all
2015-10-30T08:41:54:  Model type ................: ExtendedSource
2015-10-30T08:41:54:  Model components ..........: "EllipticalGauss" * "PowerLaw" * "Constant" 
2015-10-30T08:41:54:  Number of parameters ......: 9
2015-10-30T08:41:54:  Number of spatial par's ...: 5
2015-10-30T08:41:54:   RA .......................: 83.6201 +/- 0.00540797 [-360,360] deg (free,scale=1)
2015-10-30T08:41:54:   DEC ......................: 22.0229 +/- 0.00495831 [-90,90] deg (free,scale=1)
2015-10-30T08:41:54:   PA .......................: 45 [-360,360] deg (fixed,scale=1)
2015-10-30T08:41:54:   MajorRadius ..............: 0.198383 +/- 0.00435272 [0.001,10] deg (free,scale=1)
2015-10-30T08:41:54:   MinorRadius ..............: 0.197431 +/- 0.00433051 [0.001,10] deg (free,scale=1)
2015-10-30T08:41:54:  Number of spectral par's ..: 3
2015-10-30T08:41:54:   Prefactor ................: 5.65536e-16 +/- 1.34654e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
2015-10-30T08:41:54:   Index ....................: -2.47701 +/- 0.0183248 [-0,-5]  (free,scale=-1,gradient)
2015-10-30T08:41:54:   PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
2015-10-30T08:41:54:  Number of temporal par's ..: 1
2015-10-30T08:41:54:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
2015-10-30T08:41:54: === GCTAModelIrfBackground ===
2015-10-30T08:41:54:  Name ......................: Background
2015-10-30T08:41:54:  Instruments ...............: CTA
2015-10-30T08:41:54:  Instrument scale factors ..: unity
2015-10-30T08:41:54:  Observation identifiers ...: all
2015-10-30T08:41:54:  Model type ................: "PowerLaw" * "Constant" 
2015-10-30T08:41:54:  Number of parameters ......: 4
2015-10-30T08:41:54:  Number of spectral par's ..: 3
2015-10-30T08:41:54:   Prefactor ................: 1.00555 +/- 0.0222552 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient)
2015-10-30T08:41:54:   Index ....................: 0.00812948 +/- 0.012662 [-5,5]  (free,scale=1,gradient)
2015-10-30T08:41:54:   PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
2015-10-30T08:41:54:  Number of temporal par's ..: 1
2015-10-30T08:41:54:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)

#3 Updated by Knödlseder Jürgen over 8 years ago

Fixing the position angle but having different sizes made the problem re-appear. The semiminor axis was 0.2, semimajor axis was 0.4.

#4 Updated by Knödlseder Jürgen over 8 years ago

I was able to remove the convergence problems by limiting the theta_max to 2 times the semimajor axis value. This truncates the Gaussian quite a bit, but seriously helps the convergence. Limiting to 3 times did still produce a bunch of stalls.

So this is a kind of dirty fix as it does not guarantee that it will always work, and it may bias the fitted values somewhat. I’ll check now the pull distributions for that.

2015-10-30T09:08:03: +=================================+
2015-10-30T09:08:03: | Maximum likelihood optimisation |
2015-10-30T09:08:03: +=================================+
2015-10-30T09:08:28:  >Iteration   0: -logL=108594.263, Lambda=1.0e-03
2015-10-30T09:08:54:  >Iteration   1: -logL=108573.932, Lambda=1.0e-03, delta=20.330, max(|grad|)=-593.564853 [MinorRadius:4]
2015-10-30T09:09:19:  >Iteration   2: -logL=108569.495, Lambda=1.0e-04, delta=4.437, max(|grad|)=-209.463817 [MinorRadius:4]
2015-10-30T09:09:44:  >Iteration   3: -logL=108568.824, Lambda=1.0e-05, delta=0.671, max(|grad|)=-56.621377 [MinorRadius:4]
2015-10-30T09:10:09:  >Iteration   4: -logL=108568.769, Lambda=1.0e-06, delta=0.055, max(|grad|)=-9.553020 [MinorRadius:4]
2015-10-30T09:10:35:  >Iteration   5: -logL=108568.766, Lambda=1.0e-07, delta=0.003, max(|grad|)=-2.084455 [MajorRadius:3]
...
2015-10-30T09:11:00: +=========================================+
2015-10-30T09:11:00: | Maximum likelihood optimization results |
2015-10-30T09:11:00: +=========================================+
2015-10-30T09:11:00: === GOptimizerLM ===
2015-10-30T09:11:00:  Optimized function value ..: 108568.766
2015-10-30T09:11:00:  Absolute precision ........: 0.005
2015-10-30T09:11:00:  Acceptable value decrease .: 2
2015-10-30T09:11:00:  Optimization status .......: converged
2015-10-30T09:11:00:  Number of parameters ......: 13
2015-10-30T09:11:00:  Number of free parameters .: 9
2015-10-30T09:11:00:  Number of iterations ......: 5
2015-10-30T09:11:00:  Lambda ....................: 1e-08
2015-10-30T09:11:00:  Maximum log likelihood ....: -108568.766
2015-10-30T09:11:00:  Observed events  (Nobs) ...: 15997.000
2015-10-30T09:11:00:  Predicted events (Npred) ..: 15996.994 (Nobs - Npred = 0.00569292)
2015-10-30T09:11:00: === GModels ===
2015-10-30T09:11:00:  Number of models ..........: 2
2015-10-30T09:11:00:  Number of parameters ......: 13
2015-10-30T09:11:00: === GModelSky ===
2015-10-30T09:11:00:  Name ......................: Crab
2015-10-30T09:11:00:  Instruments ...............: all
2015-10-30T09:11:00:  Instrument scale factors ..: unity
2015-10-30T09:11:00:  Observation identifiers ...: all
2015-10-30T09:11:00:  Model type ................: ExtendedSource
2015-10-30T09:11:00:  Model components ..........: "EllipticalGauss" * "PowerLaw" * "Constant" 
2015-10-30T09:11:00:  Number of parameters ......: 9
2015-10-30T09:11:00:  Number of spatial par's ...: 5
2015-10-30T09:11:00:   RA .......................: 83.6225 +/- 0.00786815 [-360,360] deg (free,scale=1)
2015-10-30T09:11:00:   DEC ......................: 22.0224 +/- 0.0076174 [-90,90] deg (free,scale=1)
2015-10-30T09:11:00:   PA .......................: 46.9073 +/- 1.57835 [-360,360] deg (free,scale=1)
2015-10-30T09:11:00:   MajorRadius ..............: 0.421999 +/- 0.00760163 [0.001,10] deg (free,scale=1)
2015-10-30T09:11:00:   MinorRadius ..............: 0.218288 +/- 0.00384325 [0.001,10] deg (free,scale=1)
2015-10-30T09:11:00:  Number of spectral par's ..: 3
2015-10-30T09:11:00:   Prefactor ................: 6.13256e-16 +/- 1.6225e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
2015-10-30T09:11:00:   Index ....................: -2.48327 +/- 0.0194069 [-0,-5]  (free,scale=-1,gradient)
2015-10-30T09:11:00:   PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
2015-10-30T09:11:00:  Number of temporal par's ..: 1
2015-10-30T09:11:00:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
2015-10-30T09:11:00: === GCTAModelIrfBackground ===
2015-10-30T09:11:00:  Name ......................: Background
2015-10-30T09:11:00:  Instruments ...............: CTA
2015-10-30T09:11:00:  Instrument scale factors ..: unity
2015-10-30T09:11:00:  Observation identifiers ...: all
2015-10-30T09:11:00:  Model type ................: "PowerLaw" * "Constant" 
2015-10-30T09:11:00:  Number of parameters ......: 4
2015-10-30T09:11:00:  Number of spectral par's ..: 3
2015-10-30T09:11:00:   Prefactor ................: 1.04264 +/- 0.0227002 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient)
2015-10-30T09:11:00:   Index ....................: 0.0229344 +/- 0.0125059 [-5,5]  (free,scale=1,gradient)
2015-10-30T09:11:00:   PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
2015-10-30T09:11:00:  Number of temporal par's ..: 1
2015-10-30T09:11:00:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)

#5 Updated by Knödlseder Jürgen over 8 years ago

Limiting to 2.5 times the semimajor axis seams also okay. Let’s take that.

#6 Updated by Knödlseder Jürgen over 8 years ago

% Done changed from 10 to 80

It looks like the problem came from using inconsistently small angle and full trigonometric computations. As I don’t know how to do everything for full trigonometry I switched to small angle approximation and things behave much better. I selected a Gaussian cutoff of 3 sigma, as 2 sigma seems to overestimate the ellipse size.

Pull distributions will follow.

#7 Updated by Knödlseder Jürgen over 8 years ago

% Done changed from 80 to 40

There are still problems with stalls.

I switch back to a cut-off value of 2.5 sigma, which seems to behave better, but there are still stalls.

I suspect that this comes from the IRF integration method which assumes spherical coordinates instead of a small angle approximation, however I cannot confirm this at this stage. I should try to use spherical coordinates consistently, but this needs more time. Unless there is a bias, it does not need to be fixed for release 1.0.

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

Although the stalls do not disappear in all cases, I tested small and large ellipses and things look overall ok for 2.5.

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

Priority changed from Immediate to Normal
% Done changed from 40 to 50

But 2.5 looks biased in the size of the ellipse (which is too large), I therefore switched back to 3.0.

#10 Updated by Knödlseder Jürgen over 7 years ago

Target version set to 1.2.0

#11 Updated by Knödlseder Jürgen about 7 years ago

Target version changed from 1.2.0 to 1.3.0

#12 Updated by Knödlseder Jürgen almost 7 years ago

Target version changed from 1.3.0 to 1.4.0

#13 Updated by Knödlseder Jürgen over 6 years ago

Target version deleted (~~1.4.0~~)

Also available in: Atom PDF

GammaLib

Issues

Tags