Thierry Louge / gammalib

Commit fb7130403d6024bd5eed26d0c88867a0ed595e3f

Authored by Jürgen Knödlseder
1 parent 4bd3d7d4

Correct log-likelihood curvature for CTA On/Off fitting. 

Also document the formulae in the method description.
Inline Side-by-side
Showing 4 changed files with 180 additions and 50 deletions
doc/Doxyfile
  View file @fb71304
... ... @@ -256,7 +256,7 @@ MARKDOWN_SUPPORT = YES
256 256 # func(std::string) {}). This also makes the inheritance and collaboration
257 257 # diagrams that involve STL classes more complete and accurate.
258 258  
259   -BUILTIN_STL_SUPPORT = NO
  259 +BUILTIN_STL_SUPPORT = YES
260 260  
261 261 # If you use Microsoft's C++/CLI language, you should set this option to YES to
262 262 # enable parsing support.
... ... @@ -629,7 +629,7 @@ WARN_IF_DOC_ERROR = YES
629 629 # wrong or incomplete parameter documentation, but not about the absence of
630 630 # documentation.
631 631  
632   -WARN_NO_PARAMDOC = NO
  632 +WARN_NO_PARAMDOC = YES
633 633  
634 634 # The WARN_FORMAT tag determines the format of the warning messages that
635 635 # doxygen can produce. The string should contain the $file, $line, and $text
... ...
inst/cta/src/GCTAOnOffObservation.cpp
  View file @fb71304
... ... @@ -75,8 +75,7 @@ const double minerr = 1.0e-100; //!< Minimum statistical error
75 75 /* __ Coding definitions _________________________________________________ */
76 76  
77 77 /* __ Debug definitions __________________________________________________ */
78   -//#define G_EVAL_TIMING
79   -//#define G_OPT_DEBUG
  78 +//#define G_LIKELIHOOD_DEBUG //!< Debug likelihood computation
80 79  
81 80  
82 81 /*==========================================================================
... ... @@ -601,7 +600,7 @@ void GCTAOnOffObservation::compute_response(const GCTAObservation&amp; obs,
601 600 }
602 601  
603 602  
604   -/***********************************************************************
  603 +/***********************************************************************//**
605 604 * @brief Evaluate log-likelihood function for On/Off analysis
606 605 *
607 606 * @param[in] models Models.
... ... @@ -615,14 +614,92 @@ void GCTAOnOffObservation::compute_response(const GCTAObservation&amp; obs,
615 614 *
616 615 * Computes the log-likelihood value for the On/Off observation. The
617 616 * method loops over the energy bins to update the function value, its
618   - * derivatives and the curvature matrix. The number of On and Off counts
619   - * are taken from the On and Off spectra, the expected number of
620   - * gamma-ray and background events are computed from the spectral models
621   - * of the relevant components in the model container (spatial and
622   - * temporal components are ignored so far). See the model_on() and
623   - * model_off() for details about the model computations.
624   - *
625   - * @todo Add formula for curvature matrix computation
  617 + * derivatives and the curvature matrix. The number of On counts
  618 + * \f$N_{\rm on}\f$ and Off counts \f$N_{\rm off}\f$ are taken from the
  619 + * On and Off spectra, the expected number of gamma-ray events
  620 + * \f$N_{\gamma}\f$ and background events \f$N_{\rm bgd}\f$ are
  621 + * computed from the spectral models of the relevant components in the
  622 + * model container (spatial and temporal components are ignored so far).
  623 + * See the model_on() and model_off() methods for details about the
  624 + * model computations.
  625 + *
  626 + * The log-likelihood is given by
  627 + *
  628 + * \f[
  629 + * \ln L = \sum_i \ln L_i
  630 + * \f]
  631 + *
  632 + * where the sum is taken over all energy bins \f$i\f$ and
  633 + *
  634 + * \f[
  635 + * \ln L_i = - N_{\rm on} \ln N_{\rm pred} + N_{\rm pred}
  636 + * - N_{\rm off} \ln N_{\rm bgd} + N_{\rm bgd}
  637 + * \f]
  638 + *
  639 + * with
  640 + *
  641 + * \f[
  642 + * N_{\rm pred} = N_{\gamma} + \alpha N_{\rm bgd}
  643 + * \f]
  644 + *
  645 + * being the total number of predicted events for an energy bin in the On
  646 + * region,
  647 + * \f$N_{\rm on}\f$ is the total number of observed events for an energy
  648 + * bin in the On region,
  649 + * \f$N_{\rm off}\f$ is the total number of observed events for an energy
  650 + * bin in the Off region, and
  651 + * \f$N_{\rm bgd}\f$ is the predicted number of background events for an
  652 + * energy bin in the Off region.
  653 + *
  654 + * The log-likelihood gradient with respect to sky model parameters
  655 + * \f$p_{\rm sky}\f$ is given by
  656 + *
  657 + * \f[
  658 + * \left( \frac{\partial \ln L_i}{\partial p_{\rm sky}} \right) =
  659 + * \left( 1 - \frac{N_{\rm on}}{N_{\rm pred}} \right)
  660 + * \left( \frac{\partial N_{\gamma}}{\partial p_{\rm sky}} \right)
  661 + * \f]
  662 + *
  663 + * and with respect to background model parameters \f$p_{\rm bgd}\f$ is
  664 + * given by
  665 + *
  666 + * \f[
  667 + * \left( \frac{\partial \ln L_i}{\partial p_{\rm bgd}} \right) =
  668 + * \left( 1 + \alpha - \frac{N_{\rm off}}{N_{\rm bgd}} -
  669 + * \frac{\alpha N_{\rm on}}{N_{\rm pred}} \right)
  670 + * \left( \frac{\partial N_{\rm bgd}}{\partial p_{\rm bgd}} \right)
  671 + * \f]
  672 + *
  673 + * The curvature matrix elements are given by
  674 + *
  675 + * \f[
  676 + * \left( \frac{\partial^2 \ln L_i}{\partial^2 p_{\rm sky}} \right) =
  677 + * \left( \frac{N_{\rm on}}{N_{\rm pred}^2} \right)
  678 + * \left( \frac{\partial N_{\gamma}}{\partial p_{\rm sky}} \right)^2
  679 + * \f]
  680 + *
  681 + * \f[
  682 + * \left( \frac{\partial^2 \ln L_i}{\partial p_{\rm sky}
  683 + * \partial p_{\rm bgd}} \right) =
  684 + * \left( \frac{\alpha N_{\rm on}}{N_{\rm pred}^2} \right)
  685 + * \left( \frac{\partial N_{\gamma}}{\partial p_{\rm sky}} \right)
  686 + * \left( \frac{\partial N_{\rm bgd}}{\partial p_{\rm bgd}} \right)
  687 + * \f]
  688 + *
  689 + * \f[
  690 + * \left( \frac{\partial^2 \ln L_i}{\partial p_{\rm bgd}
  691 + * \partial p_{\rm sky}} \right) =
  692 + * \left( \frac{\alpha N_{\rm on}}{N_{\rm pred}^2} \right)
  693 + * \left( \frac{\partial N_{\rm bgd}}{\partial p_{\rm bgd}} \right)
  694 + * \left( \frac{\partial N_{\gamma}}{\partial p_{\rm sky}} \right)
  695 + * \f]
  696 + *
  697 + * \f[
  698 + * \left( \frac{\partial^2 \ln L_i}{\partial^2 p_{\rm bgd}} \right) =
  699 + * \left( \frac{N_{\rm off}}{N_{\rm bgd}^2} +
  700 + * \frac{\alpha^2 N_{\rm on}}{N_{\rm pred}^2} \right)
  701 + * \left( \frac{\partial N_{\rm bgd}}{\partial p_{\rm bgd}} \right)^2
  702 + * \f]
626 703 ***********************************************************************/
627 704 double GCTAOnOffObservation::likelihood(const GModels& models,
628 705 GVector* gradient,
... ... @@ -630,12 +707,12 @@ double GCTAOnOffObservation::likelihood(const GModels&amp; models,
630 707 double* npred) const
631 708 {
632 709 // Timing measurement
633   - #if defined(G_EVAL_TIMING)
  710 + #if defined(G_LIKELIHOOD_DEBUG)
634 711 clock_t t_start = clock();
635 712 #endif
636 713  
637 714 // Initialise statistics
638   - #if defined(G_OPT_DEBUG)
  715 + #if defined(G_LIKELIHOOD_DEBUG)
639 716 int n_bins = m_on_spec.size();
640 717 int n_used = 0;
641 718 int n_small_model = 0;
... ... @@ -682,7 +759,7 @@ double GCTAOnOffObservation::likelihood(const GModels&amp; models,
682 759 // Skip bin if model is too small (avoids -Inf or NaN gradients)
683 760 double nonpred = ngam + m_alpha[i] * nbgd;
684 761 if ((nbgd <= minmod) || (nonpred <= minmod)) {
685   - #if defined(G_OPT_DEBUG)
  762 + #if defined(G_LIKELIHOOD_DEBUG)
686 763 n_small_model++;
687 764 #endif
688 765 continue;
... ... @@ -695,7 +772,7 @@ double GCTAOnOffObservation::likelihood(const GModels&amp; models,
695 772 *npred += nonpred;
696 773  
697 774 // Update statistics
698   - #if defined(G_OPT_DEBUG)
  775 + #if defined(G_LIKELIHOOD_DEBUG)
699 776 n_used++;
700 777 sum_data += non;
701 778 sum_model += nonpred;
... ... @@ -712,7 +789,9 @@ double GCTAOnOffObservation::likelihood(const GModels&amp; models,
712 789 // Loop over all parameters
713 790 for (int j = 0; j < npars; ++j) {
714 791  
715   - // If spectral model for sky component is non-zero and non-infinite
  792 + // If spectral model for sky component is non-zero and
  793 + // non-infinite then handle sky component gradients and
  794 + // second derivatives including at least a sky component ...
716 795 if (sky_grad[j] != 0.0 && !gammalib::is_infinite(sky_grad[j])) {
717 796  
718 797 // Gradient
... ... @@ -722,15 +801,19 @@ double GCTAOnOffObservation::likelihood(const GModels&amp; models,
722 801 for (int k = 0; k < npars; ++k) {
723 802  
724 803 // If spectral model for sky component is non-zero and
725   - // non-infinite ...
726   - if (sky_grad[k] != 0.0 && !gammalib::is_infinite(sky_grad[k])) {
  804 + // non-infinite then we have the curvature element
  805 + // of a sky component
  806 + if (sky_grad[k] != 0.0 &&
  807 + !gammalib::is_infinite(sky_grad[k])) {
727 808 colvar[k] = sky_grad[j] * sky_grad[k] * fb;
728 809 }
729 810  
730   - // ... else if spectral model for sky component is
731   - // non-zero and non-infinite
732   - // ???????????????????? sky_grad[j]=0 ?????????????????
733   - else if (bgd_grad[k] != 0.0 && !gammalib::is_infinite(bgd_grad[k])) {
  811 + // ... else if spectral model for background component
  812 + // is non-zero and non-infinite then we have the mixed
  813 + // curvature element between a sky and a background
  814 + // component
  815 + else if (bgd_grad[k] != 0.0 &&
  816 + !gammalib::is_infinite(bgd_grad[k])) {
734 817 colvar[k] = sky_grad[j] * bgd_grad[k] * fc;
735 818 }
736 819  
... ... @@ -746,9 +829,12 @@ double GCTAOnOffObservation::likelihood(const GModels&amp; models,
746 829  
747 830 } // endif: spectral model is non-zero and non-infinite
748 831  
749   - // ... else if spectral model for background component is
750   - // non-zero and non-infinite
751   - else if (bgd_grad[j] != 0.0 && !gammalib::is_infinite(bgd_grad[j])) {
  832 + // ... otherwise if spectral model for background component is
  833 + // non-zero and non-infinite then handle background component
  834 + // gradients and second derivatives including at least a
  835 + // background component
  836 + else if (bgd_grad[j] != 0.0 &&
  837 + !gammalib::is_infinite(bgd_grad[j])) {
752 838  
753 839 // Gradient
754 840 (*gradient)[j] += bgd_factor * bgd_grad[j];
... ... @@ -756,19 +842,24 @@ double GCTAOnOffObservation::likelihood(const GModels&amp; models,
756 842 // Hessian (from first-order derivatives only)
757 843 for (int k = 0; k < npars; ++k) {
758 844  
759   - // If spectral model for sky component is non-zero and non-infinite
760   - // ???????????????????? sky_grad[j] ?????????????????
761   - if (sky_grad[k] != 0.0 && !gammalib::is_infinite(sky_grad[k])) {
762   - colvar[k] = bgd_grad[j] * bgd_grad[k] * fc;
  845 + // If spectral model for sky component is non-zero and
  846 + // non-infinite then we have the mixed curvature element
  847 + // between a sky and a background component
  848 + if (sky_grad[k] != 0.0 &&
  849 + !gammalib::is_infinite(sky_grad[k])) {
  850 + colvar[k] = bgd_grad[j] * sky_grad[k] * fc;
763 851 }
764 852  
765   - // ... else if spectral model for sky component is
766   - // non-zero and non-infinite
767   - else if (bgd_grad[k] != 0.0 && !gammalib::is_infinite(bgd_grad[k])) {
  853 + // ... else if spectral model for background component
  854 + // is non-zero and non-infinite then we have the
  855 + // curvature element of a background component
  856 + else if (bgd_grad[k] != 0.0 &&
  857 + !gammalib::is_infinite(bgd_grad[k])) {
768 858 colvar[k] = bgd_grad[j] * bgd_grad[k] * fd;
  859 + }
769 860  
770 861 // ... else neither sky nor background
771   - } else {
  862 + else {
772 863 colvar[k] = 0.0;
773 864 }
774 865  
... ... @@ -777,9 +868,8 @@ double GCTAOnOffObservation::likelihood(const GModels&amp; models,
777 868 // Update matrix
778 869 curvature->add_to_column(j, colvar);
779 870  
780   - } // endif: spectral model for background component is non-
781   - // zero and non-infinite
782   -
  871 + } // endif: spectral model for background component valid
  872 +
783 873 } // endfor: looped over all parameters for derivatives computation
784 874  
785 875 } // endfor: looped over energy bins
... ... @@ -795,7 +885,7 @@ double GCTAOnOffObservation::likelihood(const GModels&amp; models,
795 885 }
796 886  
797 887 // Dump statistics
798   - #if defined(G_OPT_DEBUG)
  888 + #if defined(G_LIKELIHOOD_DEBUG)
799 889 std::cout << "Number of parameters: " << npars << std::endl;
800 890 std::cout << "Number of bins: " << n_bins << std::endl;
801 891 std::cout << "Number of bins used for computation: " << n_used << std::endl;
... ... @@ -808,13 +898,13 @@ double GCTAOnOffObservation::likelihood(const GModels&amp; models,
808 898 #endif
809 899  
810 900 // Optionally dump gradient and curvature matrix
811   - #if defined(G_EVAL_DEBUG)
  901 + #if defined(G_LIKELIHOOD_DEBUG)
812 902 std::cout << *gradient << std::endl;
813 903 std::cout << *curvature << std::endl;
814 904 #endif
815 905  
816 906 // Timing measurement
817   - #if defined(G_EVAL_TIMING)
  907 + #if defined(G_LIKELIHOOD_DEBUG)
818 908 #ifdef _OPENMP
819 909 double t_elapse = omp_get_wtime()-t_start;
820 910 #else
... ... @@ -1488,7 +1578,7 @@ double GCTAOnOffObservation::model_off(const GModels&amp; models,
1488 1578 }
1489 1579  
1490 1580 // Fall through if model does not apply to specific instrument
1491   - // and observation identifier
  1581 + // and observation identifier.
1492 1582 if (!mptr->is_valid(this->instrument(), this->id())) {
1493 1583 ipar += mptr->size();
1494 1584 continue;
... ...
inst/cta/test/data/onoff_model.xml 0 → 100644
  View file @fb71304
  1 +<?xml version="1.0" standalone="no"?>
  2 +<source_library title="source library">
  3 + <source name="Crab" type="PointSource">
  4 + <spectrum type="PowerLaw">
  5 + <parameter name="Prefactor" scale="1e-16" value="5.7" min="1e-07" max="1000.0" free="1"/>
  6 + <parameter name="Index" scale="-1" value="2.48" min="0.0" max="+5.0" free="1"/>
  7 + <parameter name="Scale" scale="1e6" value="0.3" min="0.01" max="1000.0" free="0"/>
  8 + </spectrum>
  9 + <spatialModel type="SkyDirFunction">
  10 + <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="0"/>
  11 + <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="0"/>
  12 + </spatialModel>
  13 + </source>
  14 + <source name="Background" type="CTAIrfBackground" instrument="CTAOnOff">
  15 + <spectrum type="PowerLaw">
  16 + <parameter name="Prefactor" scale="1.0" value="1.0" min="1e-3" max="1e3" free="1"/>
  17 + <parameter name="Index" scale="1.0" value="0.0" min="-5.0" max="+5.0" free="1"/>
  18 + <parameter name="Scale" scale="1e6" value="1.0" min="0.01" max="1000.0" free="0"/>
  19 + </spectrum>
  20 + </source>
  21 +</source_library>
... ...
inst/cta/test/test_CTA.py
  View file @fb71304
... ... @@ -67,7 +67,7 @@ class Test(gammalib.GPythonTestSuite):
67 67 self.append(self.test_psf, 'Test CTA PSF classes')
68 68 self.append(self.test_edisp, 'Test CTA energy dispersion classes')
69 69 self.append(self.test_response, 'Test CTA response classes')
70   - self.append(self.test_onoff, 'Test CTA ON/OFF analysis')
  70 + self.append(self.test_onoff, 'Test CTA On/Off analysis')
71 71  
72 72 # Return
73 73 return
... ... @@ -265,7 +265,6 @@ class Test(gammalib.GPythonTestSuite):
265 265 ereco = gammalib.GEbounds(10, gammalib.GEnergy(0.1, 'TeV'),
266 266 gammalib.GEnergy(10.0, 'TeV'))
267 267  
268   -
269 268 # Create On/Off observations by filling all events found in
270 269 # the observation container and computing the response
271 270 filename = self.caldb + '/../data/irf_unbinned.xml'
... ... @@ -278,7 +277,7 @@ class Test(gammalib.GPythonTestSuite):
278 277 outobs.append(onoff)
279 278  
280 279 # Load model container and attach it to the observations
281   - models = gammalib.GModels(self.caldb + '/../data/crab_irf.xml')
  280 + models = gammalib.GModels(self.caldb + '/../data/onoff_model.xml')
282 281 outobs.models(models)
283 282  
284 283 # Perform maximum likelihood fit
... ... @@ -289,15 +288,35 @@ class Test(gammalib.GPythonTestSuite):
289 288 #print(outobs)
290 289 #print(outobs.models())
291 290  
  291 + # Test On/Off model fitting results
  292 + sky = outobs.models()['Crab']
  293 + bgd = outobs.models()['Background']
  294 + self.test_value(sky['Prefactor'].value(), 6.15281e-16, 1.0e-20,
  295 + 'Check sky model prefactor value')
  296 + self.test_value(sky['Prefactor'].error(), 2.35727e-17, 1.0e-20,
  297 + 'Check sky model prefactor error')
  298 + self.test_value(sky['Index'].value(), -2.51599, 1.0e-4,
  299 + 'Check sky model index value')
  300 + self.test_value(sky['Index'].error(), 0.037628, 1.0e-4,
  301 + 'Check sky model index error')
  302 + self.test_value(bgd['Prefactor'].value(), 1.04038, 1.0e-4,
  303 + 'Check background model prefactor value')
  304 + self.test_value(bgd['Prefactor'].error(), 0.130852, 1.0e-4,
  305 + 'Check background model prefactor error')
  306 + self.test_value(bgd['Index'].value(), 0.539376, 1.0e-4,
  307 + 'Check background model index value')
  308 + self.test_value(bgd['Index'].error(), 0.0855112, 1.0e-4,
  309 + 'Check background model index error')
  310 +
292 311 # Save PHA, ARF and RMFs
293 312 for run in outobs:
294   - run.on_spec().save('onoff_pha_on.fits', True)
295   - run.off_spec().save('onoff_pha_off.fits', True)
296   - run.arf().save('onoff_arf.fits', True)
297   - run.rmf().save('onoff_rmf.fits', True)
  313 + run.on_spec().save('test_cta_onoff_pha_on.fits', True)
  314 + run.off_spec().save('test_cta_onoff_pha_off.fits', True)
  315 + run.arf().save('test_cta_onoff_arf.fits', True)
  316 + run.rmf().save('test_cta_onoff_rmf.fits', True)
298 317  
299 318 # Save On/Off observations
300   - outobs.save('onoff.xml')
  319 + outobs.save('test_cta_onoff.xml')
301 320  
302 321 # Return
303 322 return
... ...