ROOT logo

From $ROOTSYS/tutorials/roostats/JeffreysPriorDemo.C

// tutorial demonstrating and validates the RooJeffreysPrior class

/*
JeffreysPriorDemo.C

author Kyle Cranmer
date   Dec. 2010

This tutorial demonstraites and validates the RooJeffreysPrior class

Jeffreys's prior is an 'objective prior' based on formal rules.  
It is calculated from the Fisher information matrix.  

Read more:
http://en.wikipedia.org/wiki/Jeffreys_prior

The analytic form is not known for most PDFs, but it is for 
simple cases like the Poisson mean, Gaussian mean, Gaussian sigma.

This class uses numerical tricks to calculate the Fisher Information Matrix
efficiently.  In particular, it takes advantage of a property of the
'Asimov data' as described in 
Asymptotic formulae for likelihood-based tests of new physics
Glen Cowan, Kyle Cranmer, Eilam Gross, Ofer Vitells
http://arxiv.org/abs/arXiv:1007.1727

This Demo has four parts:
TestJeffreysPriorDemo -- validates Poisson mean case 1/sqrt(mu)
TestJeffreysGaussMean -- validates Gaussian mean case
TestJeffreysGaussSigma -- validates Gaussian sigma case 1/sigma
TestJeffreysGaussMeanAndSigma -- demonstraites 2-d example

*/

#include "RooJeffreysPrior.h"

#include "RooWorkspace.h"
#include "RooDataHist.h"
#include "RooGenericPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooFitResult.h"
#include "TMatrixDSym.h"
#include "RooRealVar.h"
#include "RooAbsPdf.h"
#include "RooNumIntConfig.h"
#include "TH1F.h"

using namespace RooFit;

void JeffreysPriorDemo(){
  RooWorkspace w("w");
  w.factory("Uniform::u(x[0,1])");
  w.factory("mu[100,1,200]");
  w.factory("ExtendPdf::p(u,mu)");

  //  w.factory("Poisson::pois(n[0,inf],mu)");

  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());
  //  RooDataHist* asimov2 = w.pdf("pois")->generateBinned(*w.var("n"),ExpectedData());

  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));

  asimov->Print();
  res->Print();
  TMatrixDSym cov = res->covarianceMatrix();
  cout << "variance = " << (cov.Determinant()) << endl;
  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
  cov.Invert();
  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;

  w.defineSet("poi","mu");
  w.defineSet("obs","x");
  //  w.defineSet("obs2","n");

  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
  //  pi.specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D")  ;
  //  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",10);

  //  JeffreysPrior pi2("jeffreys2","jeffreys",*w.pdf("pois"),*w.set("poi"),*w.set("obs2"));

  //  return;
  RooGenericPdf* test = new RooGenericPdf("test","test","1./sqrt(mu)",*w.set("poi"));
  
  TCanvas* c1 = new TCanvas;
  RooPlot* plot = w.var("mu")->frame();
  //  pi.plotOn(plot, Normalization(1,RooAbsReal::Raw),Precision(.1));
  pi.plotOn(plot);
  //  pi2.plotOn(plot,LineColor(kGreen),LineStyle(kDotted));
  test->plotOn(plot,LineColor(kRed));
  plot->Draw();
 
}


//_________________________________________________
void TestJeffreysGaussMean(){
  RooWorkspace w("w");
  w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,0,10])");
  w.factory("n[10,.1,200]");
  w.factory("ExtendPdf::p(g,n)");
  w.var("sigma")->setConstant();
  w.var("n")->setConstant();

  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());

  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));

  asimov->Print();
  res->Print();
  TMatrixDSym cov = res->covarianceMatrix();
  cout << "variance = " << (cov.Determinant()) << endl;
  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
  cov.Invert();
  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;

  //  w.defineSet("poi","mu,sigma");
  w.defineSet("poi","mu");
  w.defineSet("obs","x");

  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
  //  pi.specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D")  ;
  //  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",3);

  const RooArgSet* temp = w.set("poi");
  pi.getParameters(*temp)->Print();

  //  return;
  RooGenericPdf* test = new RooGenericPdf("test","test","1",*w.set("poi"));
  
  TCanvas* c1 = new TCanvas;
  RooPlot* plot = w.var("mu")->frame();
  pi.plotOn(plot);
  test->plotOn(plot,LineColor(kRed),LineStyle(kDotted));
  plot->Draw();

  
}

//_________________________________________________
void TestJeffreysGaussSigma(){
  // this one is VERY sensitive
  // if the Gaussian is narrow ~ range(x)/nbins(x) then the peak isn't resolved
  //   and you get really bizzare shapes
  // if the Gaussian is too wide range(x) ~ sigma then PDF gets renormalized
  //   and the PDF falls off too fast at high sigma
  RooWorkspace w("w");
  w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,1,5])");
  w.factory("n[100,.1,2000]");
  w.factory("ExtendPdf::p(g,n)");
  //  w.var("sigma")->setConstant();
  w.var("mu")->setConstant();
  w.var("n")->setConstant();
  w.var("x")->setBins(301);

  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());

  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));

  asimov->Print();
  res->Print();
  TMatrixDSym cov = res->covarianceMatrix();
  cout << "variance = " << (cov.Determinant()) << endl;
  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
  cov.Invert();
  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;


  //  w.defineSet("poi","mu,sigma");
  //w.defineSet("poi","mu,sigma,n");
  w.defineSet("poi","sigma");
  w.defineSet("obs","x");

  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
  //  pi.specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D")  ;
  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",3);

  const RooArgSet* temp = w.set("poi");
  pi.getParameters(*temp)->Print();
  //  return;

  //  return;
  RooGenericPdf* test = new RooGenericPdf("test","test","sqrt(2.)/sigma",*w.set("poi"));
  
  TCanvas* c1 = new TCanvas;
  RooPlot* plot = w.var("sigma")->frame();
  pi.plotOn(plot);
  test->plotOn(plot,LineColor(kRed),LineStyle(kDotted));
  plot->Draw();

  
}


//_________________________________________________
void TestJeffreysGaussMeanAndSigma(){
  // this one is VERY sensitive
  // if the Gaussian is narrow ~ range(x)/nbins(x) then the peak isn't resolved
  //   and you get really bizzare shapes
  // if the Gaussian is too wide range(x) ~ sigma then PDF gets renormalized
  //   and the PDF falls off too fast at high sigma
  RooWorkspace w("w");
  w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,1,5])");
  w.factory("n[100,.1,2000]");
  w.factory("ExtendPdf::p(g,n)");
  //  w.var("sigma")->setConstant();
  //  w.var("mu")->setConstant();
  w.var("n")->setConstant();
  w.var("x")->setBins(301);

  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());

  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));

  asimov->Print();
  res->Print();
  TMatrixDSym cov = res->covarianceMatrix();
  cout << "variance = " << (cov.Determinant()) << endl;
  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
  cov.Invert();
  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;


  w.defineSet("poi","mu,sigma");
  //w.defineSet("poi","mu,sigma,n");
  //  w.defineSet("poi","sigma");
  w.defineSet("obs","x");

  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
  //  pi.specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D")  ;
  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",3);

  const RooArgSet* temp = w.set("poi");
  pi.getParameters(*temp)->Print();
  //  return;

  TCanvas* c1 = new TCanvas;
  TH1* Jeff2d = pi.createHistogram("2dJeffreys",*w.var("mu"),Binning(10),YVar(*w.var("sigma"),Binning(10)));
  Jeff2d->Draw("surf");
  
}


 JeffreysPriorDemo.C:1
 JeffreysPriorDemo.C:2
 JeffreysPriorDemo.C:3
 JeffreysPriorDemo.C:4
 JeffreysPriorDemo.C:5
 JeffreysPriorDemo.C:6
 JeffreysPriorDemo.C:7
 JeffreysPriorDemo.C:8
 JeffreysPriorDemo.C:9
 JeffreysPriorDemo.C:10
 JeffreysPriorDemo.C:11
 JeffreysPriorDemo.C:12
 JeffreysPriorDemo.C:13
 JeffreysPriorDemo.C:14
 JeffreysPriorDemo.C:15
 JeffreysPriorDemo.C:16
 JeffreysPriorDemo.C:17
 JeffreysPriorDemo.C:18
 JeffreysPriorDemo.C:19
 JeffreysPriorDemo.C:20
 JeffreysPriorDemo.C:21
 JeffreysPriorDemo.C:22
 JeffreysPriorDemo.C:23
 JeffreysPriorDemo.C:24
 JeffreysPriorDemo.C:25
 JeffreysPriorDemo.C:26
 JeffreysPriorDemo.C:27
 JeffreysPriorDemo.C:28
 JeffreysPriorDemo.C:29
 JeffreysPriorDemo.C:30
 JeffreysPriorDemo.C:31
 JeffreysPriorDemo.C:32
 JeffreysPriorDemo.C:33
 JeffreysPriorDemo.C:34
 JeffreysPriorDemo.C:35
 JeffreysPriorDemo.C:36
 JeffreysPriorDemo.C:37
 JeffreysPriorDemo.C:38
 JeffreysPriorDemo.C:39
 JeffreysPriorDemo.C:40
 JeffreysPriorDemo.C:41
 JeffreysPriorDemo.C:42
 JeffreysPriorDemo.C:43
 JeffreysPriorDemo.C:44
 JeffreysPriorDemo.C:45
 JeffreysPriorDemo.C:46
 JeffreysPriorDemo.C:47
 JeffreysPriorDemo.C:48
 JeffreysPriorDemo.C:49
 JeffreysPriorDemo.C:50
 JeffreysPriorDemo.C:51
 JeffreysPriorDemo.C:52
 JeffreysPriorDemo.C:53
 JeffreysPriorDemo.C:54
 JeffreysPriorDemo.C:55
 JeffreysPriorDemo.C:56
 JeffreysPriorDemo.C:57
 JeffreysPriorDemo.C:58
 JeffreysPriorDemo.C:59
 JeffreysPriorDemo.C:60
 JeffreysPriorDemo.C:61
 JeffreysPriorDemo.C:62
 JeffreysPriorDemo.C:63
 JeffreysPriorDemo.C:64
 JeffreysPriorDemo.C:65
 JeffreysPriorDemo.C:66
 JeffreysPriorDemo.C:67
 JeffreysPriorDemo.C:68
 JeffreysPriorDemo.C:69
 JeffreysPriorDemo.C:70
 JeffreysPriorDemo.C:71
 JeffreysPriorDemo.C:72
 JeffreysPriorDemo.C:73
 JeffreysPriorDemo.C:74
 JeffreysPriorDemo.C:75
 JeffreysPriorDemo.C:76
 JeffreysPriorDemo.C:77
 JeffreysPriorDemo.C:78
 JeffreysPriorDemo.C:79
 JeffreysPriorDemo.C:80
 JeffreysPriorDemo.C:81
 JeffreysPriorDemo.C:82
 JeffreysPriorDemo.C:83
 JeffreysPriorDemo.C:84
 JeffreysPriorDemo.C:85
 JeffreysPriorDemo.C:86
 JeffreysPriorDemo.C:87
 JeffreysPriorDemo.C:88
 JeffreysPriorDemo.C:89
 JeffreysPriorDemo.C:90
 JeffreysPriorDemo.C:91
 JeffreysPriorDemo.C:92
 JeffreysPriorDemo.C:93
 JeffreysPriorDemo.C:94
 JeffreysPriorDemo.C:95
 JeffreysPriorDemo.C:96
 JeffreysPriorDemo.C:97
 JeffreysPriorDemo.C:98
 JeffreysPriorDemo.C:99
 JeffreysPriorDemo.C:100
 JeffreysPriorDemo.C:101
 JeffreysPriorDemo.C:102
 JeffreysPriorDemo.C:103
 JeffreysPriorDemo.C:104
 JeffreysPriorDemo.C:105
 JeffreysPriorDemo.C:106
 JeffreysPriorDemo.C:107
 JeffreysPriorDemo.C:108
 JeffreysPriorDemo.C:109
 JeffreysPriorDemo.C:110
 JeffreysPriorDemo.C:111
 JeffreysPriorDemo.C:112
 JeffreysPriorDemo.C:113
 JeffreysPriorDemo.C:114
 JeffreysPriorDemo.C:115
 JeffreysPriorDemo.C:116
 JeffreysPriorDemo.C:117
 JeffreysPriorDemo.C:118
 JeffreysPriorDemo.C:119
 JeffreysPriorDemo.C:120
 JeffreysPriorDemo.C:121
 JeffreysPriorDemo.C:122
 JeffreysPriorDemo.C:123
 JeffreysPriorDemo.C:124
 JeffreysPriorDemo.C:125
 JeffreysPriorDemo.C:126
 JeffreysPriorDemo.C:127
 JeffreysPriorDemo.C:128
 JeffreysPriorDemo.C:129
 JeffreysPriorDemo.C:130
 JeffreysPriorDemo.C:131
 JeffreysPriorDemo.C:132
 JeffreysPriorDemo.C:133
 JeffreysPriorDemo.C:134
 JeffreysPriorDemo.C:135
 JeffreysPriorDemo.C:136
 JeffreysPriorDemo.C:137
 JeffreysPriorDemo.C:138
 JeffreysPriorDemo.C:139
 JeffreysPriorDemo.C:140
 JeffreysPriorDemo.C:141
 JeffreysPriorDemo.C:142
 JeffreysPriorDemo.C:143
 JeffreysPriorDemo.C:144
 JeffreysPriorDemo.C:145
 JeffreysPriorDemo.C:146
 JeffreysPriorDemo.C:147
 JeffreysPriorDemo.C:148
 JeffreysPriorDemo.C:149
 JeffreysPriorDemo.C:150
 JeffreysPriorDemo.C:151
 JeffreysPriorDemo.C:152
 JeffreysPriorDemo.C:153
 JeffreysPriorDemo.C:154
 JeffreysPriorDemo.C:155
 JeffreysPriorDemo.C:156
 JeffreysPriorDemo.C:157
 JeffreysPriorDemo.C:158
 JeffreysPriorDemo.C:159
 JeffreysPriorDemo.C:160
 JeffreysPriorDemo.C:161
 JeffreysPriorDemo.C:162
 JeffreysPriorDemo.C:163
 JeffreysPriorDemo.C:164
 JeffreysPriorDemo.C:165
 JeffreysPriorDemo.C:166
 JeffreysPriorDemo.C:167
 JeffreysPriorDemo.C:168
 JeffreysPriorDemo.C:169
 JeffreysPriorDemo.C:170
 JeffreysPriorDemo.C:171
 JeffreysPriorDemo.C:172
 JeffreysPriorDemo.C:173
 JeffreysPriorDemo.C:174
 JeffreysPriorDemo.C:175
 JeffreysPriorDemo.C:176
 JeffreysPriorDemo.C:177
 JeffreysPriorDemo.C:178
 JeffreysPriorDemo.C:179
 JeffreysPriorDemo.C:180
 JeffreysPriorDemo.C:181
 JeffreysPriorDemo.C:182
 JeffreysPriorDemo.C:183
 JeffreysPriorDemo.C:184
 JeffreysPriorDemo.C:185
 JeffreysPriorDemo.C:186
 JeffreysPriorDemo.C:187
 JeffreysPriorDemo.C:188
 JeffreysPriorDemo.C:189
 JeffreysPriorDemo.C:190
 JeffreysPriorDemo.C:191
 JeffreysPriorDemo.C:192
 JeffreysPriorDemo.C:193
 JeffreysPriorDemo.C:194
 JeffreysPriorDemo.C:195
 JeffreysPriorDemo.C:196
 JeffreysPriorDemo.C:197
 JeffreysPriorDemo.C:198
 JeffreysPriorDemo.C:199
 JeffreysPriorDemo.C:200
 JeffreysPriorDemo.C:201
 JeffreysPriorDemo.C:202
 JeffreysPriorDemo.C:203
 JeffreysPriorDemo.C:204
 JeffreysPriorDemo.C:205
 JeffreysPriorDemo.C:206
 JeffreysPriorDemo.C:207
 JeffreysPriorDemo.C:208
 JeffreysPriorDemo.C:209
 JeffreysPriorDemo.C:210
 JeffreysPriorDemo.C:211
 JeffreysPriorDemo.C:212
 JeffreysPriorDemo.C:213
 JeffreysPriorDemo.C:214
 JeffreysPriorDemo.C:215
 JeffreysPriorDemo.C:216
 JeffreysPriorDemo.C:217
 JeffreysPriorDemo.C:218
 JeffreysPriorDemo.C:219
 JeffreysPriorDemo.C:220
 JeffreysPriorDemo.C:221
 JeffreysPriorDemo.C:222
 JeffreysPriorDemo.C:223
 JeffreysPriorDemo.C:224
 JeffreysPriorDemo.C:225
 JeffreysPriorDemo.C:226
 JeffreysPriorDemo.C:227
 JeffreysPriorDemo.C:228
 JeffreysPriorDemo.C:229
 JeffreysPriorDemo.C:230
 JeffreysPriorDemo.C:231
 JeffreysPriorDemo.C:232
 JeffreysPriorDemo.C:233
 JeffreysPriorDemo.C:234
 JeffreysPriorDemo.C:235
 JeffreysPriorDemo.C:236
 JeffreysPriorDemo.C:237
 JeffreysPriorDemo.C:238
 JeffreysPriorDemo.C:239
 JeffreysPriorDemo.C:240
 JeffreysPriorDemo.C:241
 JeffreysPriorDemo.C:242
 JeffreysPriorDemo.C:243
 JeffreysPriorDemo.C:244