Simon Bonnefoy / ctools

Commit 43fc2cfc297d5fb6af3819c41c372bff14cfa854

Authored by Jürgen Knödlseder
1 parent 9538a7c2

Fit covariance computation in ctbutterfly (#2121)

Inline Side-by-side
Showing 4 changed files with 26 additions and 12 deletions
ChangeLog
  View file @43fc2cf
1   -2017-05-23
  1 +2017-05-31
2 2  
3 3 * Version 1.3.0 released
4 4 ========================
5 5  
  6 + Correct covariance matrix computation in ctbutterfly (#2121)
6 7 Add missing *.i SWIG file into release tarball (#2114)
7 8 ctobssim adds Monte Carlo identifiers to simulated events (#1483)
8 9 Add "phase" parameter to ctselect (#1996)
... ...
1 1 New Features and Important Changes in ctools 1.3.0.dev1
2 2  
3   -23 May 2017
  3 +31 May 2017
4 4  
5 5  
6 6 Introduction
... ... @@ -134,6 +134,9 @@ is still available for power law spectral models as an option. The "method"
134 134 parameter has been added to allow the selection of the computation method
135 135 (#1862).
136 136  
  137 +A bug in the covariance matrix has been corrected that probably was introduced
  138 +in the implementation of the Gaussian error propagation method (#2121).
  139 +
137 140  
138 141 ctulimit - Compute upper limit
139 142 ------------------------------
... ...
doc/source/users/tutorials/1dc/first_spectrum.png
View file @43fc2cf

44 KB | W: | H:

55.9 KB | W: | H:

  • 2-up
  • Swipe
  • Onion skin
src/ctbutterfly/ctbutterfly.cpp
  View file @43fc2cf
... ... @@ -243,13 +243,15 @@ void ctbutterfly::run(void)
243 243 m_obs.optimize(m_opt);
244 244 m_obs.errors(m_opt);
245 245  
246   - // Get covariance matrix.
  246 + // Get covariance matrix. For the envelope method we need an unscaled
  247 + // matrix, hence we directly access the curvature matrix and invert
  248 + // it
247 249 GObservations::likelihood likelihood = m_obs.function();
248 250 if (m_method == "GAUSSIAN") {
249   - m_covariance = likelihood.curvature()->invert();
  251 + m_covariance = likelihood.covariance();
250 252 }
251 253 else {
252   - m_covariance = likelihood.covariance();
  254 + m_covariance = likelihood.curvature()->invert();
253 255 }
254 256  
255 257 // Write optimizer, optimised model and covariance matrix into logger
... ... @@ -274,12 +276,14 @@ void ctbutterfly::run(void)
274 276 GObservations::likelihood likelihood = m_obs.function();
275 277 likelihood.eval(pars);
276 278  
277   - // Get covariance matrix.
  279 + // Get covariance matrix. For the envelope method we need an unscaled
  280 + // matrix, hence we directly access the curvature matrix and invert
  281 + // it
278 282 if (m_method == "GAUSSIAN") {
279   - m_covariance = likelihood.curvature()->invert();
  283 + m_covariance = likelihood.covariance();
280 284 }
281 285 else {
282   - m_covariance = likelihood.covariance();
  286 + m_covariance = likelihood.curvature()->invert();
283 287 }
284 288  
285 289 // Write covariance matrix into logger
... ... @@ -503,6 +507,8 @@ void ctbutterfly::get_parameters(void)
503 507  
504 508 /***********************************************************************//**
505 509 * @brief Compute butterfly using the ellipsoid boundary method
  510 + *
  511 + * @param[in] models Models.
506 512 ***************************************************************************/
507 513 void ctbutterfly::ellipsoid_boundary(GModels& models)
508 514 {
... ... @@ -560,9 +566,9 @@ void ctbutterfly::ellipsoid_boundary(GModels& models)
560 566 log_value(NORMAL, "Eigenvector 1", vector1.print());
561 567 log_value(NORMAL, "Eigenvector 2", vector2.print());
562 568  
563   - // Confidence scaling, calculated from a Chi-Sqr.
564   - // distribution with 2 dof
565   - double scale = std::sqrt( -2. * std::log(1.-m_confidence) );
  569 + // Confidence scaling, calculated from a Chi-Squared distribution
  570 + // with 2 degrees of freedom
  571 + double scale = std::sqrt(-2.0 * std::log(1.0-m_confidence));
566 572  
567 573 // Write confidence level information into logger
568 574 log_value(NORMAL, "Confidence level", m_confidence);
... ... @@ -657,6 +663,8 @@ void ctbutterfly::ellipsoid_boundary(GModels& models)
657 663 /***********************************************************************//**
658 664 * @brief Compute butterfly using Gaussian error propagation
659 665 *
  666 + * @param[in] models Models.
  667 + *
660 668 * Computes the butterfly diagram using Gaussian error propagation. This
661 669 * works for any kind of model.
662 670 ***************************************************************************/
... ... @@ -720,7 +728,9 @@ void ctbutterfly::gaussian_error_propagation(GModels& models)
720 728 // Loop over all model parameters
721 729 for (int k = 0; k < skymodel->size(); ++k) {
722 730  
723   - // Set gradient if we deal with spectral component
  731 + // Set gradient if we deal with spectral component.
  732 + // Remember that the gradient() method returns full
  733 + // parameter gradients (including the scale factor)
724 734 bool parIsSpectral = false;
725 735 for (int l = 0; l < spectral->size(); ++l) {
726 736 if ((&((*spectral)[l])) == (&((*skymodel)[k]))) {
... ...