NAME
SOOT::Examples::Math - SOOT Examples for Math
DESCRIPTION
This is a listing of all SOOT examples for Math.
EXAMPLES
FeldmanCousins.pl
use strict;
use warnings;
use SOOT ':all';
# FIXME: this SEGV's somehow...
sub _FeldmanCousins {
# Example macro of using the TFeldmanCousins class in root.
# Author : Adrian John Bevan <bevan@SLAC.Stanford.EDU>
# get a FeldmanCousins calculation object with the default limits
# of calculating a 90% CL with the minimum signal value scanned
# = 0.0 and the maximum signal value scanned of 50.0
my $f = TFeldmanCousins->new(0.9, "");
# calculate either the upper or lower limit for 10 observerd
# events with an estimated background of 3. The calculation of
# either upper or lower limit will return that limit and fill
# data members with both the upper and lower limit for you.
my $Nobserved = 10.0;
my $Nbackground = 3.0;
my $ul = $f->CalculateUpperLimit($Nobserved, $Nbackground);
my $ll = $f->GetLowerLimit();
print <<VERBATIM;
For $Nobserved data observed with and estimated background
of $Nbackground candidates, the Feldman-Cousins method of
calculating confidence limits gives:
Upper Limit = $ul
Limit = $ll
at the 90% CL
VERBATIM
}
_FeldmanCousins();binomial.pl
use strict;
use warnings;
use SOOT ':all';
sub _binomialSimple {
#
# Simple test for the binomial distribution
#
printf("\nTMath:::Binomial simple test\n");
printf("Build the Tartaglia triangle\n");
printf("============================\n");
use constant max => 13;
foreach my $i (0..max-1) {
printf "n=%2d", $i;
print " " x (max-$i);
for my $j (0..$i) {
my $bin = TMath::Nint( TMath::Binomial($i,$j));
printf("%4d", $bin);
}
print "\n";
}
}
sub _binomialFancy {
my $serr = 0;
use constant nmax => 10000;
print <<'VERBATIM';
TMath:::Binomial fancy test
Verify Newton formula for (x+y)^n
x,y in [-2,2] and n from 0 to 9
=================================
VERBATIM
my $val = 0.;
my ($x, $y);
for (0..nmax-1) {
do {
$x = 2 * (1 - 2*rand());
$y = 2 * (1 - 2*rand());
$val = abs($x+$y)*1.;
} while ($val < 0.75); # Avoid large cancellations
foreach my $j (0..9) {
my $res1 = ($x+$y) ** $j;
my $res2 = 0;
foreach my $k (0..$j) {
$res2 += $x**$k
* $y**($j-$k)
* TMath::Binomial($j,$k);
}
my $err = abs($res1-$res2)/abs($res1);
print "res1=$res1 res2=$res2 x=$x y=$y err=$err j=$j\n" if $err > 1e-10;
$serr += $err;
}
}
print "Average Error = ". $serr/nmax . "\n";
}
_binomialSimple;
_binomialFancy;limit.pl
use strict;
use warnings;
use SOOT ':all';
my $c1 = TCanvas->new("c1","Dynamic Filling Example",200,10,700,500);
$c1->SetFillColor(42);
# Create some histograms
my $background = TH1D->new("background","The expected background",30,-4,4);
my $signal = TH1D->new("signal","the expected signal",30,-4,4);
my $data = TH1D->new("data","some fake data points",30,-4,4);
$background->SetFillColor(48);
$signal->SetFillColor(41);
$data->SetMarkerStyle(21);
$data->SetMarkerColor(kBlue);
$background->Sumw2; # needed for stat uncertainty
$signal->Sumw2; # needed for stat uncertainty
# Fill histograms randomly
my $r = TRandom->new;
my ($bg, $sig, $dt);
for (0..24999) {
$bg = $r->Gaus(0.,1.)*1.0;
$sig = $r->Gaus(1.,.2)*1.0;
$background->Fill($bg,0.02);
$signal->Fill($sig,0.001);
}
for (0..499) {
$dt = $r->Gaus(0.,1.)*1.0;
$data->Fill($dt,1.0);
}
my $hs = THStack->new("hs","Signal and background compared to data...");
$hs->Add($background);
$hs->Add($signal);
$hs->Draw("hist");
$data->Draw("PE1,Same");
$c1->Modified;
$c1->Update;
my $frame = $c1->GetFrame;
$frame->SetFillColor(21);
$frame->SetBorderSize(6);
$frame->SetBorderMode(-1);
$c1->Modified;
$c1->Update;
$gSystem->ProcessEvents;
# Compute the limits
my $ds = TLimitDataSource->new($signal, $background, $data);
my $l = TLimit->new();
my $cl = $l->ComputeLimit($ds, 50000);
printCL($cl, "Computing limits...");
# Add stat uncertainty
my $scl = $l->ComputeLimit($ds,50000,1);
printCL($scl, "Computing limits with stat systematics...");
# Add some systematics
my $errorb = TH1D->new("errorb","errors on background",1,0,1);
my $errors = TH1D->new("errors","errors on signal",1,0,1);
my $names = TObjArray->new;
my $name1 = TObjString->new("bg uncertainty");
my $name2 = TObjString->new("sig uncertainty");
$names->AddLast($name1);
$names->AddLast($name2);
$errorb->SetBinContent(0,0.05); # error source 1: 5%
$errorb->SetBinContent(1,0); # error source 2: 0%
$errors->SetBinContent(0,0); # error source 1: 0%
$errors->SetBinContent(1,0.01); # error source 2: 1%
my $nds = TLimitDataSource->new;
$nds->AddChannel(
$signal, $background, $data,
TVectorD->new($errors->GetNbinsX(), $errors->GetArray()), # FIXME AddChannel expects a TVectorD argument, but that's really TVectorT<double>, which is templated and not really supported by SOOT...
TVectorD->new($errorb->GetNbinsX(), $errorb->GetArray()),
$names
);
my $ncl = $l->ComputeLimit($nds,50000,1);
printCL($ncl, "Computing limits with systematics...");
# show canonical -2lnQ plots in a new canvas
# - The histogram of -2lnQ for background hypothesis (full)
# - The histogram of -2lnQ for signal and background hypothesis (dashed)
my $c2 = TCanvas->new("c2");
$cl->Draw;
$gApplication->Run;
sub printCL {
my ($obj, $anot) = @_;
print "== ", $anot, " ==\n";
print "CLs : ", $obj->CLs, "\n";
print "CLsb : ", $obj->CLsb, "\n";
print "CLb : ", $obj->CLb, "\n";
print "<CLs> : ", $obj->GetExpectedCLs_b, "\n";
print "<CLsb> : ", $obj->GetExpectedCLsb_b, "\n";
print "<CLb> : ", $obj->GetExpectedCLb_b, "\n\n";
}mathBeta.pl
use strict;
use warnings;
use SOOT ':all';
my $c1 = TCanvas->new("c1", "TMath::BetaDist",600,800);
$c1->Divide(1, 2);
my $pad1 = $c1->cd(1);
$pad1->SetGrid();
my $fbeta = TF1->new("fbeta", "TMath::BetaDist(x, [0], [1])", 0, 1);
$fbeta->SetParameters(0.5, 0.5);
my $f1 = $fbeta->DrawCopy();
$f1->SetLineColor(kRed);
$f1->SetLineWidth(1);
$fbeta->SetParameters(0.5, 2);
my $f2 = $fbeta->DrawCopy("same");
$f2->SetLineColor(kGreen);
$f2->SetLineWidth(1);
$fbeta->SetParameters(2, 0.5);
my $f3 = $fbeta->DrawCopy("same");
$f3->SetLineColor(kBlue);
$f3->SetLineWidth(1);
$fbeta->SetParameters(2, 2);
my $f4 = $fbeta->DrawCopy("same");
$f4->SetLineColor(kMagenta);
$f4->SetLineWidth(1);
my $legend1 = TLegend->new(.5,.7,.8,.9);
$legend1->AddEntry($f1,"p=0.5 q=0.5","l");
$legend1->AddEntry($f2,"p=0.5 q=2","l");
$legend1->AddEntry($f3,"p=2 q=0.5","l");
$legend1->AddEntry($f4,"p=2 q=2","l");
$legend1->Draw();
my $pad2 = $c1->cd(2);
$pad2->SetGrid();
my $fbetai = TF1->new("fbetai", "TMath::BetaDistI(x, [0], [1])", 0, 1);
$fbetai->SetParameters(0.5, 0.5);
my $g1=$fbetai->DrawCopy();
$g1->SetLineColor(kRed);
$g1->SetLineWidth(1);
$fbetai->SetParameters(0.5, 2);
my $g2 = $fbetai->DrawCopy("same");
$g2->SetLineColor(kGreen);
$g2->SetLineWidth(1);
$fbetai->SetParameters(2, 0.5);
my $g3 = $fbetai->DrawCopy("same");
$g3->SetLineColor(kBlue);
$g3->SetLineWidth(1);
$fbetai->SetParameters(2, 2);
my $g4 = $fbetai->DrawCopy("same");
$g4->SetLineColor(kMagenta);
$g4->SetLineWidth(1);
my $legend2 = TLegend->new(.7,.15,0.9,.35);
$legend2->AddEntry($f1,"p=0.5 q=0.5","l");
$legend2->AddEntry($f2,"p=0.5 q=2","l");
$legend2->AddEntry($f3,"p=2 q=0.5","l");
$legend2->AddEntry($f4,"p=2 q=2","l");
$legend2->Draw();
$c1->cd();
$c1->Update();
$gApplication->Run;mathGammaNormal.pl
use strict;
use warnings;
use SOOT ':all';
my $myc = TCanvas->new("c1","gamma and lognormal",10,10,600,800);
$myc->Divide(1,2);
my $pad1 = $myc->cd(1);
$pad1->SetLogy();
$pad1->SetGrid();
# TMath::GammaDist
my $fgamma = TF1->new("fgamma", "TMath::GammaDist(x, [0], [1], [2])", 0, 10);
$fgamma->SetParameters(0.5, 0, 1);
my $f1 = $fgamma->DrawCopy();
$f1->SetMinimum(1e-5);
$f1->SetLineColor(kRed);
$fgamma->SetParameters(1, 0, 1);
my $f2 = $fgamma->DrawCopy("same");
$f2->SetLineColor(kGreen);
$fgamma->SetParameters(2, 0, 1);
my $f3 = $fgamma->DrawCopy("same");
$f3->SetLineColor(kBlue);
$fgamma->SetParameters(5, 0, 1);
my $f4 = $fgamma->DrawCopy("same");
$f4->SetLineColor(kMagenta);
my $legend1 = TLegend->new(.2,.15,.5,.4);
$legend1->AddEntry($f1,"gamma = 0.5 mu = 0 beta = 1","l");
$legend1->AddEntry($f2,"gamma = 1 mu = 0 beta = 1","l");
$legend1->AddEntry($f3,"gamma = 2 mu = 0 beta = 1","l");
$legend1->AddEntry($f4,"gamma = 5 mu = 0 beta = 1","l");
$legend1->Draw();
# TMath::LogNormal
my $pad2 = $myc->cd(2);
$pad2->SetLogy();
$pad2->SetGrid();
my $flog = TF1->new("flog", "TMath::LogNormal(x, [0], [1], [2])", 0, 5);
$flog->SetParameters(0.5, 0, 1);
my $g1 = $flog->DrawCopy();
$g1->SetLineColor(kRed);
$flog->SetParameters(1, 0, 1);
my $g2 = $flog->DrawCopy("same");
$g2->SetLineColor(kGreen);
$flog->SetParameters(2, 0, 1);
my $g3 = $flog->DrawCopy("same");
$g3->SetLineColor(kBlue);
$flog->SetParameters(5, 0, 1);
my $g4 = $flog->DrawCopy("same");
$g4->SetLineColor(kMagenta);
my $legend2 = TLegend->new(.2,.15,.5,.4);
$legend2->AddEntry($g1,"sigma = 0.5 theta = 0 m = 1","l");
$legend2->AddEntry($g2,"sigma = 1 theta = 0 m = 1","l");
$legend2->AddEntry($g3,"sigma = 2 theta = 0 m = 1","l");
$legend2->AddEntry($g4,"sigma = 5 theta = 0 m = 1","l");
$legend2->Draw();
$myc->Update();
$gApplication->Run;mathLaplace.pl
use strict;
use warnings;
use SOOT ':all';
# Test the TMath::LaplaceDist and TMath::LaplaceDistI functions
# author: Anna Kreshuk
my $c1 = TCanvas->new("c1", "TMath::LaplaceDist",600,800);
$c1->Divide(1, 2);
my $pad1 = $c1->cd(1);
$pad1->SetGrid();
my $flaplace = TF1->new("flaplace", "TMath::LaplaceDist(x, [0], [1])", -10, 10);
$flaplace->SetParameters(0, 1);
my $f1 = $flaplace->DrawCopy();
$f1->SetLineColor(kRed);
$f1->SetLineWidth(1);
$flaplace->SetParameters(0, 2);
my $f2 = $flaplace->DrawCopy("same");
$f2->SetLineColor(kGreen);
$f2->SetLineWidth(1);
$flaplace->SetParameters(2, 1);
my $f3 = $flaplace->DrawCopy("same");
$f3->SetLineColor(kBlue);
$f3->SetLineWidth(1);
$flaplace->SetParameters(2, 2);
my $f4 = $flaplace->DrawCopy("same");
$f4->SetLineColor(kMagenta);
$f4->SetLineWidth(1);
my $legend1 = TLegend->new(.7,.7,.9,.9);
$legend1->AddEntry($f1,"alpha=0 beta=1","l");
$legend1->AddEntry($f2,"alpha=0 beta=2","l");
$legend1->AddEntry($f3,"alpha=2 beta=1","l");
$legend1->AddEntry($f4,"alpha=2 beta=2","l");
$legend1->Draw();
my $pad2 = $c1->cd(2);
$pad2->SetGrid();
my $flaplacei = TF1->new("flaplacei", "TMath::LaplaceDistI(x, [0], [1])", -10, 10);
$flaplacei->SetParameters(0, 1);
my $g1 = $flaplacei->DrawCopy();
$g1->SetLineColor(kRed);
$g1->SetLineWidth(1);
$flaplacei->SetParameters(0, 2);
my $g2 = $flaplacei->DrawCopy("same");
$g2->SetLineColor(kGreen);
$g2->SetLineWidth(1);
$flaplacei->SetParameters(2, 1);
my $g3 = $flaplacei->DrawCopy("same");
$g3->SetLineColor(kBlue);
$g3->SetLineWidth(1);
$flaplacei->SetParameters(2, 2);
my $g4 = $flaplacei->DrawCopy("same");
$g4->SetLineColor(kMagenta);
$g4->SetLineWidth(1);
my $legend2 = TLegend->new(.7,.15,0.9,.35);
$legend2->AddEntry($f1,"alpha=0 beta=1","l");
$legend2->AddEntry($f2,"alpha=0 beta=2","l");
$legend2->AddEntry($f3,"alpha=2 beta=1","l");
$legend2->AddEntry($f4,"alpha=2 beta=2","l");
$legend2->Draw();
$c1->cd();
$gApplication->Run;mathStudent.pl
use strict;
use warnings;
use SOOT ':all';
my $DistCanvas = TCanvas->new("DistCanvas", "Distribution graphs", 10, 10, 1000, 800);
$DistCanvas->SetFillColor(17);
$DistCanvas->Divide(2, 2);
my $pad1 = $DistCanvas->cd(1);
$pad1->SetGrid();
$pad1->SetFrameFillColor(19);
my $leg = TLegend->new(0.6, 0.7, 0.89, 0.89);
my $fgaus = TF1->new("gaus", "TMath::Gaus(x, [0], [1], [2])", -5, 5);
$fgaus->SetTitle("Student density");
$fgaus->SetLineStyle(2);
$fgaus->SetLineWidth(1);
$fgaus->SetParameters(0, 1, 1);
$leg->AddEntry($fgaus->DrawCopy(), "Normal(0,1)", "l");
my $student = TF1->new("student", "TMath::Student(x,[0])", -5, 5);
$student->SetTitle("Student density");
$student->SetLineWidth(1);
$student->SetParameter(0, 10);
$student->SetLineColor(4);
$leg->AddEntry($student->DrawCopy("lsame"), "10 degrees of freedom", "l");
$student->SetParameter(0, 3);
$student->SetLineColor(2);
$leg->AddEntry($student->DrawCopy("lsame"), "3 degrees of freedom", "l");
$student->SetParameter(0, 1);
$student->SetLineColor(1);
$leg->AddEntry($student->DrawCopy("lsame"), "1 degree of freedom", "l");
$leg->Draw();
#drawing the set of student cumulative probability functions
my $pad2 = $DistCanvas->cd(2);
$pad2->SetFrameFillColor(19);
$pad2->SetGrid();
my $studentI = TF1->new("studentI", "TMath::StudentI(x, [0])", -5, 5);
$studentI->SetTitle("Student cumulative dist.");
$studentI->SetLineWidth(1);
my $leg2 = TLegend->new(0.6, 0.4, 0.89, 0.6);
$studentI->SetParameter(0, 10);
$studentI->SetLineColor(4);
$leg2->AddEntry($studentI->DrawCopy(), "10 degrees of freedom", "l");
$studentI->SetParameter(0, 3);
$studentI->SetLineColor(2);
$leg2->AddEntry($studentI->DrawCopy("lsame"), "3 degrees of freedom", "l");
$studentI->SetParameter(0, 1);
$studentI->SetLineColor(1);
$leg2->AddEntry($studentI->DrawCopy("lsame"), "1 degree of freedom", "l");
$leg2->Draw();
#drawing the set of F-dist. densities
my $fDist = TF1->new("fDist", "TMath::FDist(x, [0], [1])", 0, 2);
$fDist->SetTitle("F-Dist. density");
$fDist->SetLineWidth(1);
my $legF1 = TLegend->new(0.7, 0.7, 0.89, 0.89);
my $pad3 = $DistCanvas->cd(3);
$pad3->SetFrameFillColor(19);
$pad3->SetGrid();
$fDist->SetParameters(1, 1);
$fDist->SetLineColor(1);
$legF1->AddEntry($fDist->DrawCopy(), "N=1 M=1", "l");
$fDist->SetParameter(1, 10);
$fDist->SetLineColor(2);
$legF1->AddEntry($fDist->DrawCopy("lsame"), "N=1 M=10", "l");
$fDist->SetParameters(10, 1);
$fDist->SetLineColor(8);
$legF1->AddEntry($fDist->DrawCopy("lsame"), "N=10 M=1", "l");
$fDist->SetParameters(10, 10);
$fDist->SetLineColor(4);
$legF1->AddEntry($fDist->DrawCopy("lsame"), "N=10 M=10", "l");
$legF1->Draw();
# drawing the set of F cumulative dist.functions
my $fDistI = TF1->new("fDist", "TMath::FDistI(x, [0], [1])", 0, 2);
$fDistI->SetTitle("Cumulative dist. function for F");
$fDistI->SetLineWidth(1);
my $legF2 = TLegend->new(0.7, 0.3, 0.89, 0.5);
my $pad4 = $DistCanvas->cd(4);
$pad4->SetFrameFillColor(19);
$pad4->SetGrid();
$fDistI->SetParameters(1, 1);
$fDistI->SetLineColor(1);
$legF2->AddEntry($fDistI->DrawCopy(), "N=1 M=1", "l");
$fDistI->SetParameters(1, 10);
$fDistI->SetLineColor(2);
$legF2->AddEntry($fDistI->DrawCopy("lsame"), "N=1 M=10", "l");
$fDistI->SetParameters(10, 1);
$fDistI->SetLineColor(8);
$legF2->AddEntry($fDistI->DrawCopy("lsame"), "N=10 M=1", "l");
$fDistI->SetParameters(10, 10);
$fDistI->SetLineColor(4);
$legF2->AddEntry($fDistI->DrawCopy("lsame"), "N=10 M=10", "l");
$legF2->Draw();
$DistCanvas->Update();
$gApplication->Run;mathcoreCDF.pl
use strict;
use SOOT ':all';
# Example macro describing how to use the different cumulative
# distribution functions in ROOT. The macro shows four of them with
# respect to their two variables. In order to run the macro type:
# FIXME doesn't work... didn't work in CINT either. Need to figure out MathCore better
$gSystem->Load("root/libMathCore");
my $f1a = TF2->new("f1a", "ROOT::Math::breitwigner_prob(x, y)", -10, 10, 0, 10);
my $f2a = TF2->new("f2a", "ROOT::Math::cauchy_quant(x,y)", 0, 20, 0,20);
my $f3a = TF2->new("f3a", "ROOT::Math::normal_quant(x,y)", -10, 10, 0, 5);
my $f4a = TF2->new("f4a", "ROOT::Math::exponential_prob(x,y)", 0, 10, 0, 5);
my $c1 = TCanvas->new("c1","c1",1000,750);
$c1->Divide(2,2);
$c1->cd(1);
$f1a->Draw("surf1");
$c1->cd(2);
$f2a->Draw("surf1");
$c1->cd(3);
$f3a->Draw("surf1");
$c1->cd(4);
$f4a->Draw("surf1");
$c1->Update;
$gApplication->Run;mathcoreSpecFunc.pl
use strict;
use warnings;
use SOOT ':all';
# Example macro describing how to use the special mathematical functions
# taking full advantage of the precision and speed of the C99 compliant
# environments. To execute the macro type in:
#
# root[0]: .x mathcoreSpecFunc.C
#
# It will create two canvases:
#
# a) one with the representation of the tgamma, lgamma, erf and erfc functions
# b) one with the relative difference between the old ROOT versions and the
# C99 implementation (on obsolete platform+compiler combinations which are
# not C99 compliant it will call the original ROOT implementations, hence
# the difference will be 0)
#
# The naming and numbering of the functions is taken from
# <A HREF="http:#www.open-std.org/jtc1/sc22/wg21/docs/papers/2004/n1687.pdf">
# Matt Austern,
# (Draft) Technical Report on Standard Library Extensions,
# N1687=04-0127, September 10, 2004</A>
#
# Author: Andras Zsenei
$gSystem->Load("libMathCore");
my $f1a = TF1->new("f1a","ROOT::Math::tgamma(x)",0,100);
my $f1b = TF1->new("f1b","TMath::Abs((ROOT::Math::tgamma(x)-TMath::Gamma(x))/ROOT::Math::tgamma(x))",0,100);
my $f2a = TF1->new("f2a","ROOT::Math::lgamma(x)",0,100);
my $f2b = TF1->new("f2b","TMath::Abs((ROOT::Math::lgamma(x)-TMath::LnGamma(x))/ROOT::Math::lgamma(x))",0,100);
my $f3a = TF1->new("f3a","ROOT::Math::erf(x)",0,5);
my $f3b = TF1->new("f3b","TMath::Abs((ROOT::Math::erf(x)-TMath::Erf(x))/ROOT::Math::erf(x))",0,5);
my $f4a = TF1->new("f4a","ROOT::Math::erfc(x)",0,5);
my $f4b = TF1->new("f4b","TMath::Abs((ROOT::Math::erfc(x)-TMath::Erfc(x))/ROOT::Math::erfc(x))",0,5);
my $c1 = TCanvas->new("c1","c1",1000,750);
$c1->Divide(2,2);
$c1->cd(1);
$f1a->Draw();
$c1->cd(2);
$f2a->Draw();
$c1->cd(3);
$f3a->Draw();
$c1->cd(4);
$f4a->Draw();
my $c2 = TCanvas->new("c2","c2",1000,750);
$c2->Divide(2,2);
$c2->cd(1);
$f1b->Draw();
$c2->cd(2);
$f2b->Draw();
$c2->cd(3);
$f3b->Draw();
$c2->cd(4);
$f4b->Draw();
$gApplication->Run;mathcoreStatFunc.pl
use strict;
use warnings;
use SOOT ':all';
# Example macro describing how to use the different probability
# density functions in ROOT. The macro shows four of them with
# respect to their two variables. In order to run the macro type:
#
# Author: Andras Zsenei
$gSystem->Load("libMathCore");
my $f1a = TF2->new("f1a","ROOT::Math::cauchy_pdf(x, y)",0,10,0,10);
my $f2a = TF2->new("f2a","ROOT::Math::chisquared_pdf(x,y)",0,20, 0,20);
my $f3a = TF2->new("f3a","ROOT::Math::gaussian_pdf(x,y)",0,10,0,5);
my $f4a = TF2->new("f4a","ROOT::Math::tdistribution_pdf(x,y)",0,10,0,5);
my $c1 = TCanvas->new("c1","c1",1000,750);
$c1->Divide(2,2);
$c1->cd(1);
$f1a->Draw("surf1");
$c1->cd(2);
$f2a->Draw("surf1");
$c1->cd(3);
$f3a->Draw("surf1");
$c1->cd(4);
$f4a->Draw("surf1");
$gApplication->Run;principal.pl
use strict;
use warnings;
use SOOT ':all';
# Principal Components Analysis (PCA) example
#
# Example of using TPrincipal as a stand alone class.
#
# We create n-dimensional data points, where c = trunc(n / 5) + 1
# are correlated with the rest n - c randomly distributed variables.
#
my $n = shift || 10;
my $m = shift || 10000;
my $c = $n / 5 + 1;
printf "*************************************************\n";
printf "* Principal Component Analysis *\n";
printf "* *\n";
printf "* Number of variables: %8d *\n", $n;
printf "* Number of data points: %8d *\n", $m;
printf "* Number of dependent variables: %4d *\n", $c;
printf "* *\n";
printf "*************************************************\n";
# Initilase the TPrincipal object. Use the empty string for the
# final argument, if you don't wan't the covariance
# matrix. Normalising the covariance matrix is a good idea if your
# variables have different orders of magnitude.
my $principal = TPrincipal->new($n,"N");
# Use a pseudo-random number generator
my $random = TRandom->new;
# Make the m data-points
# Make a variable to hold our data
# Allocate memory for the data point
my $data = [];
for my $i (0..$m-1) {
# First we create the un-correlated, random variables, according
# to one of three distributions
for my $j (0..($n-$c-1)) {
if ($j % 3 == 0) {
$data->[$j] = $random->Gaus(5,1);
}
elsif ($j % 3 == 1) {
$data->[$j] = $random->Poisson(8);
}
else {
$data->[$j] = $random->Exp(2);
}
}
# Then we create the correlated variables
for my $j (0..$c-1) {
$data->[$n - $c + $j] = 0;
for my $k (0..($n-$c-$j-1)) {
$data->[$n - $c + $j] += $data->[$k];
}
}
# Finally we're ready to add this datapoint to the PCA
$principal->AddRow($data);
}
# Do the actual analysis
$principal->MakePrincipals();
# Print out the result on
$principal->Print();
# Test the PCA
$principal->Test();
# Make some histograms of the orginal, principal, residue, etc data
$principal->MakeHistograms();
# Make two functions to map between feature and pattern space
$principal->MakeCode();
# Start a browser, so that we may browse the histograms generated
# above
my $b = TBrowser->new("principalBrowser", $principal);
$gApplication->Run;rolke.pl
use strict;
use warnings;
use SOOT ':all';
my $bm = 0.0;
my $tau = 2.5;
my $mid = 1;
my $m = 100;
my $z = 50;
my $y = 10;
my $x = 5;
# Initialize parameters not used.
my $e = 0.0;
my $em = 0.0;
my $sde=0.0;
my $sdb=0.0;
my $b = 0.0;
my $g = TRolke->new();
$g->SetCL(0.90);
my $ul = $g->CalculateInterval($x,$y,$z,$bm,$em,$e,$mid,$sde,$sdb,$tau,$b,$m);
my $ll = $g->GetLowerLimit();
print "Assuming MODEL 1\n";
print "the Profile Likelihood interval is :\n";
print "[", $ll, ",", $ul, "]\n";
$tau = 2.5;
$mid = 2;
$y = 3;
$x = 10;
$em=0.9;
$sde=0.05;
$g->SetCL(0.95);
$ul = $g->CalculateInterval($x,$y,$z,$bm,$em,$e,$mid,$sde,$sdb,$tau,$b,$m);
$ll = $g->GetLowerLimit();
print "Assuming MODEL 2\n";
print "the Profile Likelihood interval is :\n";
print "[", $ll, ",", $ul, "]\n";
$mid = 3;
$bm = 5.0;
$x = 10;
$em = 0.9;
$sde=0.05;
$sdb=0.5;
$g->SetCL(0.99);
$ul = $g->CalculateInterval($x,$y,$z,$bm,$em,$e,$mid,$sde,$sdb,$tau,$b,$m);
$ll = $g->GetLowerLimit();
print "Assuming MODEL 3\n";
print "the Profile Likelihood interval is :\n";
print "[", $ll, ",", $ul, "]\n";
$tau = 5;
$mid = 4;
$y = 7;
$x = 1;
$e = 0.25;
$g->SetCL(0.68);
$ul = $g->CalculateInterval($x,$y,$z,$bm,$em,$e,$mid,$sde,$sdb,$tau,$b,$m);
$ll = $g->GetLowerLimit();
print "Assuming MODEL 4\n";
print "the Profile Likelihood interval is :\n";
print "[", $ll, ",", $ul, "]\n";
$mid = 5;
$bm = 0.0;
$x = 1;
$e = 0.65;
$sdb=1.0;
$g->SetCL(0.80);
$ul = $g->CalculateInterval($x,$y,$z,$bm,$em,$e,$mid,$sde,$sdb,$tau,$b,$m);
$ll = $g->GetLowerLimit();
print "Assuming MODEL 5\n";
print "the Profile Likelihood interval is :\n";
print "[", $ll, ",", $ul, "]\n";
$y = 1;
$mid = 6;
$m = 750;
$z = 500;
$x = 25;
$b = 10.0;
$g->SetCL(0.90);
$ul = $g->CalculateInterval($x,$y,$z,$bm,$em,$e,$mid,$sde,$sdb,$tau,$b,$m);
$ll = $g->GetLowerLimit();
print "Assuming MODEL 6\n";
print "the Profile Likelihood interval is :\n";
print "[", $ll, ",", $ul, "]\n";
$mid = 7;
$x = 15;
$em = 0.77;
$sde=0.15;
$b = 10.0;
$g->SetCL(0.95);
$ul = $g->CalculateInterval($x,$y,$z,$bm,$em,$e,$mid,$sde,$sdb,$tau,$b,$m);
$ll = $g->GetLowerLimit();
print "Assuming MODEL 7\n";
print "the Profile Likelihood interval is :\n";
print "[", $ll, ",", $ul, "]\n";vavilov.pl
use strict;
use warnings;
use SOOT ':all';
#test of the TMath::Vavilov distribution
use constant n => 1000;
my $x = [];
my $y1 = [];
my $y2 = [];
my $y3 = [];
my $y4 = [];
my $r = TRandom->new();
for my $i (0..n-1) {
my $rv = $r->Uniform(-2, 10);
push @$x, $rv;
push @$y1, TMath::Vavilov($x->[$i], 0.30, 0.5);
push @$y2, TMath::Vavilov($x->[$i], 0.15, 0.5);
push @$y3, TMath::Vavilov($x->[$i], 0.25, 0.5);
push @$y4, TMath::Vavilov($x->[$i], 0.05, 0.5);
}
my $c1 = TCanvas->new("c1", "Vavilov density");
$c1->SetGrid();
$c1->SetHighLightColor(19);
my $gr1 = TGraph->new(n, $x, $y1);
my $gr2 = TGraph->new(n, $x, $y2);
my $gr3 = TGraph->new(n, $x, $y3);
my $gr4 = TGraph->new(n, $x, $y4);
$gr1->SetTitle("TMath::Vavilov density");
$gr1->Draw("ap");
$gr2->Draw("psame");
$gr2->SetMarkerColor(kRed);
$gr3->Draw("psame");
$gr3->SetMarkerColor(kBlue);
$gr4->Draw("psame");
$gr4->SetMarkerColor(kGreen);
$gApplication->Run;SEE ALSO
AUTHOR
Steffen Mueller, <smueller@cpan.org>
COPYRIGHT AND LICENSE
Copyright (C) 2010 by Steffen Mueller
SOOT, the Perl-ROOT wrapper, is free software; you can redistribute it and/or modify it under the same terms as ROOT itself, that is, the GNU Lesser General Public License. A copy of the full license text is available from the distribution as the LICENSE file.