#include <boost/math/distributions/inverse_chi_squared.hpp>

#include <iostream>
using namespace std;

#include <boost/math/distributions/inverse_chi_squared.hpp>
#include <boost/math/distributions/chi_squared.hpp>
//using namespace boost::math;  // Avoid this is real life!
// It will get you into trouble due to name clashes.


int main() 
{
  double events, lo, hi;

  using boost::math::inverse_chi_squared_distribution;
  using boost::math::chi_squared;
  using boost::math::quantile;

  
  events=17; 

  // Johnson original 
  inverse_chi_squared_distribution<> dist_hi(2*events);
  hi = 1 / (quantile(complement(dist_hi, 0.975)) * 2);

  inverse_chi_squared_distribution<> dist_lo(2*(events+1));
  lo = 1 / (quantile(complement(dist_lo, 0.025)) * 2);

  cout << "events = " << events << endl;
  cout << hi << endl;  // 25.983
  cout << lo << endl; //   10.6679
  //  Agreed - This program produces results of [10.67, 25.98] for 17 events.

  // John suggestion 

  //Excel  chiinv(0.025, 34) = 19.80625294 / 2 = 9.90312647
  

  lo = quantile(chi_squared(2*events),  (0.05)/2) / 2; // (34, 0.025)
  cout << "lo   " << lo << endl; // 9.90313

  hi = quantile(chi_squared(2*events), (0.95)/2) / 2; //  = (34, 0. 925);
  cout << "hi   " << hi << endl; //   hi   16.4136  
  // Above should be 1 - 0.05/2) = 0.975

    hi = quantile(chi_squared(2*events), 0.975) / 2; //  = (34, 0. 925);
  cout << "hi   " << hi << endl; //   hi   
  // = chi_squared(34, 0.025); //    hi   25.983
  // =CHISQ.INV(0.975,34)/2 = 25.9829976


     hi = quantile(chi_squared(2*(events+1)), 0.975) / 2; //  = (36, 0. 925);
  cout << "hi   " << hi << endl; //   hi   
    
   return 0;
}

/*
   chisqr_Johnson.cpp
  chisqr_eg.vcxproj -> J:\Cpp\chisqr_eg\Debug\chisqr_eg.exe
  events = 17
  25.983
  10.6679
  lo   9.90313
  hi   16.4136
  hi   25.983
  hi   27.2186


*/