test_stanfit_dummy.C
Go to the documentation of this file.
1 /*
2  * test_stanfit_dummy.C:
3  * Test the use of the MCMC tool Stan for fitting in CAFAna.
4  * This version doesn't use any neutrinos, just the fitting mechanics.
5  *
6  * Original author: J. Wolcott <jwolcott@fnal.gov>
7  * Date: November 2018
8  */
9 
10 #include <memory>
11 
12 #include "TCanvas.h"
13 #include "TGraph.h"
14 #include "TH1D.h"
15 #include "TMath.h"
16 
17 #include "CAFAna/Fit/StanConfig.h"
18 #include "CAFAna/Fit/StanFitter.h"
19 #include "CAFAna/Core/SystShifts.h"
21 #include "Utilities/func/MathUtil.h"
22 #include "Utilities/func/Stan.h"
23 
24 using namespace ana;
25 
26 namespace test
27 {
28  std::size_t N_POINTS = 5;
29  double TRUE_A = 1.8;
30 
31  // try a starting point pretty far away from the truth
32  double SEED_VAL = 25000;
33 
34  // a 'systematic' for fitting.
35  // will actually be the only free parameter in our model...
36  class QuadraticParameter : public ISyst
37  {
38  public:
40  : ISyst("QuadraticParameter", "quadratic scale factor")
41  {}
42 
43  };
45 
46  double LogGauss(double x, double mu, double sigma)
47  {
48  const double offset = log(1./sqrt(2*TMath::Pi()));
49  return offset - util::sqr(x - mu) / (2 * util::sqr(sigma));
50  }
51 
52  // a very simple experiment that uses a log-Gaussian likelihood
53  // to fit a single parameter in a quadratic function:
54  // y = a * x^2
56  {
57  public:
59  const SystShifts& syst = SystShifts::Nominal()) const override
60  {
61  stan::math::var ll = 0;
62  for (unsigned int i = 1; i <= N_POINTS; i++)
63  {
64  auto ll_contrib = stan::math::normal_lpdf(TRUE_A * util::sqr(i),
65  syst.GetShift<stan::math::var>(&kQuadParam) * util::sqr(i), 1.);
66  ll += ll_contrib;
67 // std::cout << "For a = " << syst.GetShift(&kQuadParam) << ", LL = " << ll_contrib << std::endl;
68  }
69  return ll;
70  }
71 
72  };
73 }
74 
76 {
78 
79 
80 
81  StanFitter fitter(&expt, {}, {&test::kQuadParam});
83  config.num_warmup = 1000;
84  config.num_samples = 10000;
85 // config.max_depth = 17;
86 // config.verbosity = StanConfig::Verbosity::kVerbose;
87  config.verbosity = StanConfig::Verbosity::kQuiet;
88  fitter.SetStanConfig(config);
89 
90 
91  SystShifts shifts(&test::kQuadParam, 1);
92 
93  TCanvas c;
94  TGraph ll;
95  for (int x = -25; x <= 25; x++)
96  {
97  shifts.ResetToNominal();
98  shifts.SetShift(&test::kQuadParam, x);
99 // std::cout << "x = " << x << "; chi2 = " << expt.ChiSq(nullptr, shifts)/-2 << std::endl;
100  ll.SetPoint(ll.GetN(), x, expt.LogLikelihood(nullptr, shifts).val());
101  }
102  ll.Draw("alp");
103  c.SaveAs("test_stanfit_LL.png");
104 
105  std::cout << "Seeding fit parameter at value: " << test::SEED_VAL << std::endl;
106  shifts.ResetToNominal();
108 
109  fitter.Fit(shifts);
110  //fitter.TestGradients(nullptr, shifts);
111 
112  std::cout << "Best fit: a = " << shifts.GetShift(&test::kQuadParam) << std::endl;
113  std::cout << " (true value was a = " << test::TRUE_A << ")" << std::endl;
114 
115  c.Clear();
116  const_cast<TTree*>(fitter.GetSamples().ToTTree())->Draw("QuadraticParameter");
117  c.SaveAs("test_stanfit_samples.png");
118 }
tree Draw("slc.nhit")
double TRUE_A
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
T sqrt(T number)
Definition: d0nt_math.hpp:156
Configuration parameters for the Stan MCMC fitter, bundled up together for easy storage and parameter...
Definition: StanConfig.h:6
static SystShifts Nominal()
Definition: SystShifts.h:34
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
stan::math::var LogLikelihood(osc::IOscCalcAdjustableStan *, const SystShifts &syst=SystShifts::Nominal()) const override
expt
Definition: demo5.py:34
double SEED_VAL
T GetShift(const ISyst *syst) const
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
double LogGauss(double x, double mu, double sigma)
return_type< T_y, T_loc, T_scale >::type normal_lpdf(const T_y &y, const T_loc &mu, const T_scale &sigma)
Definition: normal_lpdf.hpp:45
void test_stanfit_dummy()
void ResetToNominal()
Definition: SystShifts.cxx:143
General interface to any calculator that lets you set the parameters.
Base class defining interface for experiments.
Definition: IExperiment.h:14
double val() const
Definition: var.hpp:294
const QuadraticParameter kQuadParam
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
Fitter type that bolts the Stan fitting tools onto CAFAna.
Definition: StanFitter.h:126
std::size_t N_POINTS