Logo ROOT  
Reference Guide
RooUnbinnedL.cxx
Go to the documentation of this file.
1/*
2 * Project: RooFit
3 * Authors:
4 * PB, Patrick Bos, Netherlands eScience Center, p.bos@esciencecenter.nl
5 *
6 * Copyright (c) 2021, CERN
7 *
8 * Redistribution and use in source and binary forms,
9 * with or without modification, are permitted according to the terms
10 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)
11 */
12
13/**
14\file RooUnbinnedL.cxx
15\class RooUnbinnedL
16\ingroup Roofitcore
17
18Class RooUnbinnedL implements a -log(likelihood) calculation from a dataset
19(assumed to be unbinned) and a PDF. The NLL is calculated as
20\f[
21 \sum_\mathrm{data} -\log( \mathrm{pdf}(x_\mathrm{data}))
22\f]
23In extended mode, a
24\f$ N_\mathrm{expect} - N_\mathrm{observed}*log(N_\mathrm{expect}) \f$ term is added.
25**/
26
28
29#include "RooAbsData.h"
30#include "RooAbsPdf.h"
31#include "RooAbsDataStore.h"
32#include "RooNLLVar.h" // RooNLLVar::ComputeScalar
33#include "RunContext.h"
34#include "RooChangeTracker.h"
35
36namespace RooFit {
37namespace TestStatistics {
38
40 bool useBatchedEvaluations)
41 : RooAbsL(RooAbsL::ClonePdfData{pdf, data}, data->numEntries(), 1, extended),
42 useBatchedEvaluations_(useBatchedEvaluations)
43{
44 std::unique_ptr<RooArgSet> params(pdf->getParameters(data));
45 paramTracker_ = std::make_unique<RooChangeTracker>("chtracker","change tracker",*params,true);
46}
47
49 : RooAbsL(other), apply_weight_squared(other.apply_weight_squared), _first(other._first),
50 useBatchedEvaluations_(other.useBatchedEvaluations_), lastSection_(other.lastSection_),
51 cachedResult_(other.cachedResult_)
52{
53 paramTracker_ = std::make_unique<RooChangeTracker>(*other.paramTracker_);
54}
55
57
58//////////////////////////////////////////////////////////////////////////////////
59
60/// Returns true if value was changed, false otherwise.
62{
63 if (apply_weight_squared != flag) {
65 return true;
66 }
67 // setValueDirty();
68 return false;
69}
70
73}
74
75//////////////////////////////////////////////////////////////////////////////////
76/// Calculate and return likelihood on subset of data from firstEvent to lastEvent
77/// processed with a step size of 'stepSize'. If this an extended likelihood and
78/// and the zero event is processed the extended term is added to the return
79/// likelihood.
80///
82RooUnbinnedL::evaluatePartition(Section events, std::size_t /*components_begin*/, std::size_t /*components_end*/)
83{
84 // Throughout the calculation, we use Kahan's algorithm for summing to
85 // prevent loss of precision - this is a factor four more expensive than
86 // straight addition, but since evaluating the PDF is usually much more
87 // expensive than that, we tolerate the additional cost...
89 double sumWeight;
90
91 // Do not reevaluate likelihood if parameters nor event range have changed
92 if (!paramTracker_->hasChanged(true) && events == lastSection_ && (cachedResult_ != 0)) return cachedResult_;
93
94 data_->store()->recalculateCache(nullptr, events.begin(N_events_), events.end(N_events_), 1, true);
95
97 std::unique_ptr<RooBatchCompute::RunContext> evalData;
98 std::tie(result, sumWeight) = RooNLLVar::computeBatchedFunc(pdf_.get(), data_.get(), evalData, normSet_.get(), apply_weight_squared,
99 1, events.begin(N_events_), events.end(N_events_));
100 } else {
101 std::tie(result, sumWeight) = RooNLLVar::computeScalarFunc(pdf_.get(), data_.get(), normSet_.get(), apply_weight_squared,
102 1, events.begin(N_events_), events.end(N_events_));
103 }
104
105 // include the extended maximum likelihood term, if requested
106 if (extended_ && events.begin_fraction == 0) {
108
109 // TODO: the following should also be factored out into free/static functions like RooNLLVar::Compute*
110 // Calculate sum of weights-squared here for extended term
111 double sumW2;
113 const RooSpan<const double> eventWeights = data_->getWeightBatch(0, N_events_);
114 if (eventWeights.empty()) {
115 sumW2 = (events.end(N_events_) - events.begin(N_events_)) * data_->weightSquared();
116 } else {
118 for (std::size_t i = 0; i < eventWeights.size(); ++i) {
119 kahanWeight.AddIndexed(eventWeights[i] * eventWeights[i], i);
120 }
121 sumW2 = kahanWeight.Sum();
122 }
123 } else { // scalar mode
124 ROOT::Math::KahanSum<double> sumW2KahanSum;
125 for (Int_t i = 0; i < data_->numEntries(); i++) {
126 data_->get(i);
127 sumW2KahanSum += data_->weightSquared();
128 }
129 sumW2 = sumW2KahanSum.Sum();
130 }
131
132 double expected = pdf_->expectedEvents(data_->get());
133
134 // Adjust calculation of extended term with W^2 weighting: adjust poisson such that
135 // estimate of Nexpected stays at the same value, but has a different variance, rescale
136 // both the observed and expected count of the Poisson with a factor sum[w] / sum[w^2] which is
137 // the effective weight of the Poisson term.
138 // i.e. change Poisson(Nobs = sum[w]| Nexp ) --> Poisson( sum[w] * sum[w] / sum[w^2] | Nexp * sum[w] /
139 // sum[w^2] ) weighted by the effective weight sum[w^2]/ sum[w] in the likelihood. Since here we compute
140 // the likelihood with the weight square we need to multiply by the square of the effective weight expectedW
141 // = expected * sum[w] / sum[w^2] : effective expected entries observedW = sum[w] * sum[w] / sum[w^2] :
142 // effective observed entries The extended term for the likelihood weighted by the square of the weight will
143 // be then:
144 // (sum[w^2]/ sum[w] )^2 * expectedW - (sum[w^2]/ sum[w] )^2 * observedW * log (expectedW) and this is
145 // using the previous expressions for expectedW and observedW
146 // sum[w^2] / sum[w] * expected - sum[w^2] * log (expectedW)
147 // and since the weights are constants in the likelihood we can use log(expected) instead of log(expectedW)
148
149 double expectedW2 = expected * sumW2 / data_->sumEntries();
150 double extra = expectedW2 - sumW2 * log(expected);
151
152 // double y = pdf->extendedTerm(sumW2, data->get()) - carry;
153
154 result += extra;
155 } else {
156 result += pdf_->extendedTerm(data_->sumEntries(), data_->get());
157 }
158 }
159
160 // If part of simultaneous PDF normalize probability over
161 // number of simultaneous PDFs: -sum(log(p/n)) = -sum(log(p)) + N*log(n)
162 if (sim_count_ > 1) {
163 result += sumWeight * log(1.0 * sim_count_);
164 }
165
166 // At the end of the first full calculation, wire the caches
167 if (_first) {
168 _first = false;
169 pdf_->wireAllCaches();
170 }
171
173 lastSection_ = events;
174 return result;
175}
176
177} // namespace TestStatistics
178} // namespace RooFit
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
T Sum() const
Definition: Util.h:240
void AddIndexed(T input, std::size_t index)
Add input to the sum.
Definition: Util.h:231
RooArgSet * getParameters(const RooAbsData *data, bool stripDisconnected=true) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
Definition: RooAbsArg.cxx:541
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:62
std::shared_ptr< RooAbsData > data_
Definition: RooAbsL.h:135
std::unique_ptr< RooArgSet > normSet_
Pointer to set with observables used for normalization.
Definition: RooAbsL.h:136
std::shared_ptr< RooAbsPdf > pdf_
Definition: RooAbsL.h:134
ROOT::Math::KahanSum< double > cachedResult_
Definition: RooUnbinnedL.h:53
bool setApplyWeightSquared(bool flag)
Returns true if value was changed, false otherwise.
ROOT::Math::KahanSum< double > evaluatePartition(Section events, std::size_t components_begin, std::size_t components_end) override
Calculate and return likelihood on subset of data from firstEvent to lastEvent processed with a step ...
std::unique_ptr< RooChangeTracker > paramTracker_
Definition: RooUnbinnedL.h:51
bool apply_weight_squared
Apply weights squared?
Definition: RooUnbinnedL.h:48
RooUnbinnedL(RooAbsPdf *pdf, RooAbsData *data, RooAbsL::Extended extended=RooAbsL::Extended::Auto, bool useBatchedEvaluations=false)
static RooNLLVar::ComputeResult computeScalarFunc(const RooAbsPdf *pdfClone, RooAbsData *dataClone, RooArgSet *normSet, bool weightSq, std::size_t stepSize, std::size_t firstEvent, std::size_t lastEvent)
Definition: RooNLLVar.cxx:472
static RooNLLVar::ComputeResult computeBatchedFunc(const RooAbsPdf *pdfClone, RooAbsData *dataClone, std::unique_ptr< RooBatchCompute::RunContext > &evalData, RooArgSet *normSet, bool weightSq, std::size_t stepSize, std::size_t firstEvent, std::size_t lastEvent)
Definition: RooNLLVar.cxx:364
A simple container to hold a batch of data values.
Definition: RooSpan.h:34
constexpr std::span< T >::index_type size() const noexcept
Definition: RooSpan.h:121
constexpr bool empty() const noexcept
Definition: RooSpan.h:125
RVec< PromoteType< T > > log(const RVec< T > &v)
Definition: RVec.hxx:1748
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: Common.h:18
Convenience wrapper class used to distinguish between pdf/data owning and non-owning constructors.
Definition: RooAbsL.h:39
A part of some range delimited by two fractional points between 0 and 1 (inclusive).
Definition: RooAbsL.h:58
std::size_t begin(std::size_t N_total) const
Definition: RooAbsL.h:66
std::size_t end(std::size_t N_total) const
Definition: RooAbsL.h:68