//+ example of sampling a multi-dim distribution using the DistSampler class //NOTE: This tutorial must be run with ACLIC //Author L. Moneta Dec 2010 // function (a 4d gaussian) #include "TMath.h" #include "TF2.h" #include "TStopwatch.h" #include "Math/DistSampler.h" #include "Math/DistSamplerOptions.h" #include "Math/MinimizerOptions.h" #include "Math/Factory.h" #include "TKDTreeBinning.h" #include "TTree.h" #include "TFile.h" #include "TMatrixDSym.h" #include "TVectorD.h" #include "TCanvas.h" #include <cmath> // Gauss ND function // make a class in order to avoid constructing the // matrices for every call // This however requires that the code must be compiled with ACLIC bool debug = false; struct GausND { TVectorD X; TVectorD Mu; TMatrixDSym CovMat; GausND( int dim ) : X(TVectorD(dim)), Mu(TVectorD(dim)), CovMat(TMatrixDSym(dim) ) {} double operator() (double *x, double *p) { // 4 parameters int dim = X.GetNrows(); int k = 0; for (int i = 0; i<dim; ++i) { X[i] = x[i] - p[k]; k++; } for (int i = 0; i<dim; ++i) { CovMat(i,i) = p[k]*p[k]; k++; } for (int i = 0; i<dim; ++i) { for (int j = i+1; j<dim; ++j) { // p now are the correlations N(N-1)/2 CovMat(i,j) = p[k]*sqrt(CovMat(i,i)*CovMat(j,j)); CovMat(j,i) = CovMat(i,j); k++; } } if (debug) { X.Print(); CovMat.Print(); } double det = CovMat.Determinant(); if (det <= 0) { Fatal("GausND","Determinant is <= 0 det = %f",det); CovMat.Print(); return 0; } double norm = std::pow( 2. * TMath::Pi(), dim/2) * sqrt(det); // compute the gaussians CovMat.Invert(); double fval = std::exp( - 0.5 * CovMat.Similarity(X) )/ norm; if (debug) { std::cout << "det " << det << std::endl; std::cout << "norm " << norm << std::endl; std::cout << "fval " << fval << std::endl; } return fval; } }; using namespace ROOT::Math; void multidimSampling() { #ifdef __CINT__ std::cout << "DO NOT RUN WITH CINT:" << std::endl; std::cout << "we are using a custom function which requires" << std::endl; std::cout << "that this tutorial must be compiled with ACLIC" << std::endl; return; #endif const int N = 100000; //const int NBin = 1000; const int DIM = 4; double xmin[] = {-10,-10,-10, -10}; double xmax[] = { 10, 10, 10, 10}; double par0[] = { 1., -1., 2, 0, // the gaussian mu 1, 2, 1, 3, // the sigma 0.5,0.,0.,0.,0.,0.8 }; // the correlation const int NPAR = DIM + DIM*(DIM+1)/2; // 14 in the 4d case // generate the sample GausND gaus4d(4); TF1 * f = new TF1("functionND",gaus4d,0,1,14); f->SetParameters(par0); double x0[] = {0,0,0,0}; // for debugging if (debug) f->EvalPar(x0,0); debug = false; TString name; for (int i = 0; i < NPAR; ++i ) { if (i < DIM) f->SetParName(i, name.Format("mu_%d",i+1) ); else if (i < 2*DIM) f->SetParName(i, name.Format("sig_%d",i-DIM+1) ); else if (i < 2*DIM) f->SetParName(i, name.Format("sig_%d",i-2*DIM+1) ); } //ROOT::Math::DistSamplerOptions::SetDefaultSampler("Foam"); DistSampler * sampler = Factory::CreateDistSampler(); if (sampler == 0) { Info("multidimSampling","Default sampler %s is not available try with Foam ", ROOT::Math::DistSamplerOptions::DefaultSampler().c_str() ); ROOT::Math::DistSamplerOptions::SetDefaultSampler("Foam"); } sampler = Factory::CreateDistSampler(); if (sampler == 0) { Error("multidimSampling","Foam sampler is not available - exit "); return; } sampler->SetFunction(*f,DIM); sampler->SetRange(xmin,xmax); bool ret = sampler->Init(); std::vector<double> data1(DIM*N); double v[DIM]; TStopwatch w; if (!ret) { Error("Sampler::Init","Error initializing unuran sampler"); return; } // generate the data w.Start(); for (int i = 0; i < N; ++i) { sampler->Sample(v); for (int j = 0; j < DIM; ++j) data1[N*j + i] = v[j]; } w.Stop(); w.Print(); // fill tree with data TFile * file = new TFile("multiDimSampling.root","RECREATE"); double x[DIM]; TTree * t1 = new TTree("t1","Tree from Unuran"); t1->Branch("x",x,"x[4]/D"); for (int i = 0; i < N; ++i) { for (int j = 0; j < DIM; ++j) { x[j] = data1[i+N*j]; } t1->Fill(); } // try to fit the 2d data; // GausND gaus2(2); // TF2 * f2 = new TF2("f2",gaus2,-3,3,-3,3,5,"GausND"); // f2->SetParameters(0,0,1,1,0); // f2->SetParLimits(4,-1,1); // t1->UnbinnedFit("f2","x[0]:x[1]"); // plot the data t1->Draw("x[0]:x[1]:x[2]:x[3]","","candle"); TCanvas * c2 = new TCanvas(); c2->Divide(3,2); int ic=1; c2->cd(ic++); t1->Draw("x[0]:x[1]"); c2->cd(ic++); t1->Draw("x[0]:x[2]"); c2->cd(ic++); t1->Draw("x[0]:x[3]"); c2->cd(ic++); t1->Draw("x[1]:x[2]"); c2->cd(ic++); t1->Draw("x[1]:x[3]"); c2->cd(ic++); t1->Draw("x[2]:x[3]"); t1->Write(); file->Close(); } multidimSampling.C:1 multidimSampling.C:2 multidimSampling.C:3 multidimSampling.C:4 multidimSampling.C:5 multidimSampling.C:6 multidimSampling.C:7 multidimSampling.C:8 multidimSampling.C:9 multidimSampling.C:10 multidimSampling.C:11 multidimSampling.C:12 multidimSampling.C:13 multidimSampling.C:14 multidimSampling.C:15 multidimSampling.C:16 multidimSampling.C:17 multidimSampling.C:18 multidimSampling.C:19 multidimSampling.C:20 multidimSampling.C:21 multidimSampling.C:22 multidimSampling.C:23 multidimSampling.C:24 multidimSampling.C:25 multidimSampling.C:26 multidimSampling.C:27 multidimSampling.C:28 multidimSampling.C:29 multidimSampling.C:30 multidimSampling.C:31 multidimSampling.C:32 multidimSampling.C:33 multidimSampling.C:34 multidimSampling.C:35 multidimSampling.C:36 multidimSampling.C:37 multidimSampling.C:38 multidimSampling.C:39 multidimSampling.C:40 multidimSampling.C:41 multidimSampling.C:42 multidimSampling.C:43 multidimSampling.C:44 multidimSampling.C:45 multidimSampling.C:46 multidimSampling.C:47 multidimSampling.C:48 multidimSampling.C:49 multidimSampling.C:50 multidimSampling.C:51 multidimSampling.C:52 multidimSampling.C:53 multidimSampling.C:54 multidimSampling.C:55 multidimSampling.C:56 multidimSampling.C:57 multidimSampling.C:58 multidimSampling.C:59 multidimSampling.C:60 multidimSampling.C:61 multidimSampling.C:62 multidimSampling.C:63 multidimSampling.C:64 multidimSampling.C:65 multidimSampling.C:66 multidimSampling.C:67 multidimSampling.C:68 multidimSampling.C:69 multidimSampling.C:70 multidimSampling.C:71 multidimSampling.C:72 multidimSampling.C:73 multidimSampling.C:74 multidimSampling.C:75 multidimSampling.C:76 multidimSampling.C:77 multidimSampling.C:78 multidimSampling.C:79 multidimSampling.C:80 multidimSampling.C:81 multidimSampling.C:82 multidimSampling.C:83 multidimSampling.C:84 multidimSampling.C:85 multidimSampling.C:86 multidimSampling.C:87 multidimSampling.C:88 multidimSampling.C:89 multidimSampling.C:90 multidimSampling.C:91 multidimSampling.C:92 multidimSampling.C:93 multidimSampling.C:94 multidimSampling.C:95 multidimSampling.C:96 multidimSampling.C:97 multidimSampling.C:98 multidimSampling.C:99 multidimSampling.C:100 multidimSampling.C:101 multidimSampling.C:102 multidimSampling.C:103 multidimSampling.C:104 multidimSampling.C:105 multidimSampling.C:106 multidimSampling.C:107 multidimSampling.C:108 multidimSampling.C:109 multidimSampling.C:110 multidimSampling.C:111 multidimSampling.C:112 multidimSampling.C:113 multidimSampling.C:114 multidimSampling.C:115 multidimSampling.C:116 multidimSampling.C:117 multidimSampling.C:118 multidimSampling.C:119 multidimSampling.C:120 multidimSampling.C:121 multidimSampling.C:122 multidimSampling.C:123 multidimSampling.C:124 multidimSampling.C:125 multidimSampling.C:126 multidimSampling.C:127 multidimSampling.C:128 multidimSampling.C:129 multidimSampling.C:130 multidimSampling.C:131 multidimSampling.C:132 multidimSampling.C:133 multidimSampling.C:134 multidimSampling.C:135 multidimSampling.C:136 multidimSampling.C:137 multidimSampling.C:138 multidimSampling.C:139 multidimSampling.C:140 multidimSampling.C:141 multidimSampling.C:142 multidimSampling.C:143 multidimSampling.C:144 multidimSampling.C:145 multidimSampling.C:146 multidimSampling.C:147 multidimSampling.C:148 multidimSampling.C:149 multidimSampling.C:150 multidimSampling.C:151 multidimSampling.C:152 multidimSampling.C:153 multidimSampling.C:154 multidimSampling.C:155 multidimSampling.C:156 multidimSampling.C:157 multidimSampling.C:158 multidimSampling.C:159 multidimSampling.C:160 multidimSampling.C:161 multidimSampling.C:162 multidimSampling.C:163 multidimSampling.C:164 multidimSampling.C:165 multidimSampling.C:166 multidimSampling.C:167 multidimSampling.C:168 multidimSampling.C:169 multidimSampling.C:170 multidimSampling.C:171 multidimSampling.C:172 multidimSampling.C:173 multidimSampling.C:174 multidimSampling.C:175 multidimSampling.C:176 multidimSampling.C:177 multidimSampling.C:178 multidimSampling.C:179 multidimSampling.C:180 multidimSampling.C:181 multidimSampling.C:182 multidimSampling.C:183 multidimSampling.C:184 multidimSampling.C:185 multidimSampling.C:186 multidimSampling.C:187 multidimSampling.C:188 multidimSampling.C:189 multidimSampling.C:190 multidimSampling.C:191 multidimSampling.C:192 multidimSampling.C:193 multidimSampling.C:194 multidimSampling.C:195 multidimSampling.C:196 multidimSampling.C:197 multidimSampling.C:198 multidimSampling.C:199 multidimSampling.C:200 multidimSampling.C:201 multidimSampling.C:202 multidimSampling.C:203 multidimSampling.C:204 multidimSampling.C:205 |
|