Maximum likelihood and least squares

入门阅读：ROOT Fit 教程。这里保留少计数与高计数样本的专题比较。

Consider a decay-time measurement in which each recorded time follows an exponential distribution. The same physical parameter, the lifetime $\tau$, can be estimated from the event values themselves or from a histogram of those values. This example keeps the original sequence from unbinned likelihood to binned Poisson likelihood, least squares, and finally a high-statistics comparison.

Choose Python / PyROOT or ROOT C++ below. Both versions use the same model, random seed, fit range, and plots.

Code language: Python / PyROOT ROOT C++

Which method should be used?

Available data and model	Method	Use it when
Individual event values $t_j$ and a normalized PDF	unbinned likelihood	The event positions are available and their information should not be reduced to histogram counts. The PDF must be normalized over the analysis range.
Counts $n_i$ in histogram bins and expected counts $\mu_i$	binned Poisson likelihood	The data are counts, especially when some bins have few or zero events. In TH1::Fit, option L selects this method.
Histogram or graph values with meaningful, approximately Gaussian uncertainties	least squares	The Gaussian approximation is reasonable. For a sparse count histogram, the default $\chi^2$ fit can exclude empty bins and should not be the first choice.
Event values or bin counts together with a model for the expected total yield	extended likelihood	The observed total count carries information about a yield or rate. If the sample size is fixed by construction and only the shape matters, a non-extended likelihood is sufficient.

Binning discards the position of an event inside its bin. It is therefore a choice in the statistical analysis, not only a choice of drawing style.

Generate the same decay samples

The true lifetime is $\tau=1$. A sample with 50 events illustrates sparse histogram counts; a second sample with 10,000 events shows the high-statistics limit. Both histograms use 20 bins over $0\leq t\leq10$.

import math
from array import array
import ROOT

ROOT.gStyle.SetOptStat(0)

tau_true = 1.0
n_low = 50
n_high = 10000
rng = ROOT.TRandom3(12345)  # Fixed seed: the example is reproducible.

# Store event values for unbinned likelihood and fill histograms
# from those same events for the binned fits.
low_times = [rng.Exp(tau_true) for _ in range(n_low)]
high_times = [rng.Exp(tau_true) for _ in range(n_high)]

h_low = ROOT.TH1D("h_low_py", ";decay time t;counts / bin", 20, 0, 10)
h_high = ROOT.TH1D("h_high_py", ";decay time t;counts / bin", 20, 0, 10)
for t in low_times:
    h_low.Fill(t)
for t in high_times:
    h_high.Fill(t)

#include "TCanvas.h"
#include "TF1.h"
#include "TGraph.h"
#include "TH1D.h"
#include "TLegend.h"
#include "TLine.h"
#include "TLatex.h"
#include "TRandom3.h"
#include "TStyle.h"
#include "TFitResultPtr.h"
#include <cmath>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <vector>

gStyle->SetOptStat(0);

const double tauTrue = 1.0;
const int nLow = 50;
const int nHigh = 10000;
TRandom3 rng(12345);  // Fixed seed: the example is reproducible.

// Store event values for unbinned likelihood and fill histograms
// from those same events for the binned fits.
std::vector<double> lowTimes;
std::vector<double> highTimes;
lowTimes.reserve(nLow);
highTimes.reserve(nHigh);

auto hLow = new TH1D("hLow", ";decay time t;counts / bin", 20, 0, 10);
auto hHigh = new TH1D("hHigh", ";decay time t;counts / bin", 20, 0, 10);
for (int i = 0; i < nLow; ++i) {
    const double t = rng.Exp(tauTrue);
    lowTimes.push_back(t);
    hLow->Fill(t);
}
for (int i = 0; i < nHigh; ++i) {
    const double t = rng.Exp(tauTrue);
    highTimes.push_back(t);
    hHigh->Fill(t);
}

Unbinned likelihood

For measured times $t_1,\ldots,t_N$,

$$ L(\tau)=\prod_{j=1}^{N}\frac{1}{\tau}e^{-t_j/\tau}, \qquad \log L(\tau)=-N\log\tau-\frac{1}{\tau}\sum_j t_j . $$

The maximum occurs at the sample mean,

$$ \hat\tau=\frac{1}{N}\sum_j t_j . $$

For this model, $\mathrm{Var}(\hat\tau)=\tau^2/N$, while the standard deviation of the estimator is $\sigma_{\hat\tau}=\tau/\sqrt{N}$. Variance and standard deviation must not be confused.

The curve below shows $\Delta\log L=\log L(\tau)-\log L(\hat\tau)$. For one fitted parameter, the intersections with $\Delta\log L=-0.5$ give the usual likelihood-based one-standard-deviation interval when the quadratic approximation is adequate.

sum_low = sum(low_times)
tau_hat_low = sum_low / n_low
sigma_tau_low = tau_hat_low / math.sqrt(n_low)

tau_scan = [0.45 + i * (1.10 / 220) for i in range(221)]
logl = [-n_low * math.log(tau) - sum_low / tau for tau in tau_scan]
logl_max = -n_low * math.log(tau_hat_low) - sum_low / tau_hat_low
delta_logl = [value - logl_max for value in logl]

g_logl = ROOT.TGraph(len(tau_scan), array("d", tau_scan), array("d", delta_logl))
g_logl.SetTitle("Unbinned likelihood;#tau;#Delta log L")
g_logl.SetLineWidth(2)

c_unbinned = ROOT.TCanvas("c_unbinned_py", "Unbinned likelihood", 700, 500)
g_logl.Draw("AL")
g_logl.SetMinimum(-10)
g_logl.SetMaximum(0.5)
line_half = ROOT.TLine(tau_scan[0], -0.5, tau_scan[-1], -0.5)
line_half.SetLineStyle(2)
line_half.Draw()
line_hat = ROOT.TLine(tau_hat_low, -10, tau_hat_low, 0)
line_hat.SetLineColor(ROOT.kRed)
line_hat.Draw()
c_unbinned.Draw()

print(f"unbinned: tau = {tau_hat_low:.4f} +/- {sigma_tau_low:.4f}")

const double sumLow = std::accumulate(lowTimes.begin(), lowTimes.end(), 0.0);
const double tauHatLow = sumLow / nLow;
const double sigmaTauLow = tauHatLow / std::sqrt(nLow);

auto gLogL = new TGraph();
const int nScan = 221;
const double tauMin = 0.45;
const double tauMax = 1.55;
const double logLMax = -nLow * std::log(tauHatLow) - sumLow / tauHatLow;
for (int i = 0; i < nScan; ++i) {
    const double tau = tauMin + (tauMax - tauMin) * i / (nScan - 1);
    const double logL = -nLow * std::log(tau) - sumLow / tau;
    gLogL->SetPoint(i, tau, logL - logLMax);
}
gLogL->SetTitle("Unbinned likelihood;#tau;#Delta log L");
gLogL->SetLineWidth(2);

auto cUnbinned = new TCanvas("cUnbinned", "Unbinned likelihood", 700, 500);
gLogL->Draw("AL");
gLogL->SetMinimum(-10);
gLogL->SetMaximum(0.5);
auto lineHalf = new TLine(tauMin, -0.5, tauMax, -0.5);
lineHalf->SetLineStyle(2);
lineHalf->Draw();
auto lineHat = new TLine(tauHatLow, -10, tauHatLow, 0);
lineHat->SetLineColor(kRed);
lineHat->Draw();
cUnbinned->Draw();

std::cout << std::fixed << std::setprecision(4)
          << "unbinned: tau = " << tauHatLow
          << " +/- " << sigmaTauLow << std::endl;

unbinned: tau = 0.8799 +/- 0.1244

Binned Poisson likelihood

After the events are put into bins, the data are the counts $n_i$. If the expected count in bin $i$ is $\mu_i(\theta)$, the appropriate count model is

$$ L(\theta)=\prod_i \operatorname{Poisson}\!\left(n_i\mid\mu_i(\theta)\right). $$

The TH1::Fit options used here are:

• L: use the binned Poisson likelihood instead of the default $\chi^2$;
• I: integrate the function over each bin instead of evaluating it only at the bin center;
• R: use the range defined by the TF1;
• S: return a TFitResultPtr, which gives access to status, parameters, errors, and the covariance matrix.

Option I is useful when the function changes appreciably across a bin. It does not recover the event positions lost by binning.

h_low_poisson = h_low.Clone("h_low_poisson_py")
f_low_poisson = ROOT.TF1("f_low_poisson_py", "[0]*exp(-x/[1])", 0, 10)
f_low_poisson.SetParNames("A", "tau")
f_low_poisson.SetParameters(h_low.GetMaximum(), 1.0)
f_low_poisson.SetParLimits(0, 0, 1.0e9)
f_low_poisson.SetParLimits(1, 0.05, 5.0)
f_low_poisson.SetLineColor(ROOT.kRed)

# L: Poisson likelihood, I: bin integral, R: TF1 range,
# S: return fit result, Q: quiet, 0: do not draw automatically.
result_low_poisson = h_low_poisson.Fit(f_low_poisson, "LIRSQ0")

c_low_poisson = ROOT.TCanvas("c_low_poisson_py", "Low statistics: Poisson likelihood", 700, 500)
c_low_poisson.SetLogy()
h_low_poisson.SetMinimum(0.3)
h_low_poisson.SetTitle("50 events: binned Poisson likelihood;decay time t;counts / bin")
h_low_poisson.Draw("E")
f_low_poisson.Draw("same")
label_low_poisson = ROOT.TLatex()
label_low_poisson.SetNDC()
label_low_poisson.DrawLatex(
    0.53, 0.82,
    f"#tau = {f_low_poisson.GetParameter(1):.3f} #pm {f_low_poisson.GetParError(1):.3f}",
)
c_low_poisson.Draw()

print(
    f"binned Poisson, N=50: tau = {f_low_poisson.GetParameter(1):.4f} "
    f"+/- {f_low_poisson.GetParError(1):.4f}, "
    f"status = {int(result_low_poisson)}, ndf = {f_low_poisson.GetNDF()}"
)

auto hLowPoisson = static_cast<TH1D*>(hLow->Clone("hLowPoisson"));
auto fLowPoisson = new TF1("fLowPoisson", "[0]*exp(-x/[1])", 0, 10);
fLowPoisson->SetParNames("A", "tau");
fLowPoisson->SetParameters(hLow->GetMaximum(), 1.0);
fLowPoisson->SetParLimits(0, 0, 1.0e9);
fLowPoisson->SetParLimits(1, 0.05, 5.0);
fLowPoisson->SetLineColor(kRed);

// L: Poisson likelihood, I: bin integral, R: TF1 range,
// S: return fit result, Q: quiet, 0: do not draw automatically.
TFitResultPtr resultLowPoisson = hLowPoisson->Fit(fLowPoisson, "LIRSQ0");

auto cLowPoisson = new TCanvas("cLowPoisson", "Low statistics: Poisson likelihood", 700, 500);
cLowPoisson->SetLogy();
hLowPoisson->SetMinimum(0.3);
hLowPoisson->SetTitle("50 events: binned Poisson likelihood;decay time t;counts / bin");
hLowPoisson->Draw("E");
fLowPoisson->Draw("same");
auto labelLowPoisson = new TLatex();
labelLowPoisson->SetNDC();
labelLowPoisson->DrawLatex(
    0.53, 0.82,
    Form("#tau = %.3f #pm %.3f", fLowPoisson->GetParameter(1), fLowPoisson->GetParError(1)));
cLowPoisson->Draw();

std::cout << "binned Poisson, N=50: tau = "
          << fLowPoisson->GetParameter(1) << " +/- " << fLowPoisson->GetParError(1)
          << ", status = " << static_cast<int>(resultLowPoisson)
          << ", ndf = " << fLowPoisson->GetNDF() << std::endl;

binned Poisson, N=50: tau = 0.8970 +/- 0.1287, status = 0, ndf = 18

Least squares with sparse counts

The default histogram fit minimizes a $\chi^2$ built from the bin contents and their errors. For an unweighted count histogram, ROOT normally assigns $\sigma_i=\sqrt{n_i}$; an empty bin then has zero error and is excluded from the default $\chi^2$ fit. In the 50-event histogram, the empty tail bins contain information—the expected count should also be small—but the default least-squares fit does not use them.

This does not mean least squares is generally wrong. It is appropriate when the fitted values have known, approximately Gaussian errors. The issue here is using that approximation for sparse Poisson counts.

h_low_ls = h_low.Clone("h_low_ls_py")
f_low_ls = ROOT.TF1("f_low_ls_py", "[0]*exp(-x/[1])", 0, 10)
f_low_ls.SetParNames("A", "tau")
f_low_ls.SetParameters(h_low.GetMaximum(), 1.0)
f_low_ls.SetParLimits(0, 0, 1.0e9)
f_low_ls.SetParLimits(1, 0.05, 5.0)
f_low_ls.SetLineColor(ROOT.kBlue + 1)

# Without L, TH1::Fit uses the default chi-square fit.
result_low_ls = h_low_ls.Fit(f_low_ls, "IRSQ0")

c_low_ls = ROOT.TCanvas("c_low_ls_py", "Low statistics: least squares", 700, 500)
c_low_ls.SetLogy()
h_low_ls.SetMinimum(0.3)
h_low_ls.SetTitle("50 events: least squares;decay time t;counts / bin")
h_low_ls.Draw("E")
f_low_ls.Draw("same")
label_low_ls = ROOT.TLatex()
label_low_ls.SetNDC()
label_low_ls.DrawLatex(
    0.53, 0.82,
    f"#tau = {f_low_ls.GetParameter(1):.3f} #pm {f_low_ls.GetParError(1):.3f}",
)
c_low_ls.Draw()

print(
    f"least squares, N=50: tau = {f_low_ls.GetParameter(1):.4f} "
    f"+/- {f_low_ls.GetParError(1):.4f}, "
    f"status = {int(result_low_ls)}, ndf = {f_low_ls.GetNDF()}"
)

auto hLowLS = static_cast<TH1D*>(hLow->Clone("hLowLS"));
auto fLowLS = new TF1("fLowLS", "[0]*exp(-x/[1])", 0, 10);
fLowLS->SetParNames("A", "tau");
fLowLS->SetParameters(hLow->GetMaximum(), 1.0);
fLowLS->SetParLimits(0, 0, 1.0e9);
fLowLS->SetParLimits(1, 0.05, 5.0);
fLowLS->SetLineColor(kBlue + 1);

// Without L, TH1::Fit uses the default chi-square fit.
TFitResultPtr resultLowLS = hLowLS->Fit(fLowLS, "IRSQ0");

auto cLowLS = new TCanvas("cLowLS", "Low statistics: least squares", 700, 500);
cLowLS->SetLogy();
hLowLS->SetMinimum(0.3);
hLowLS->SetTitle("50 events: least squares;decay time t;counts / bin");
hLowLS->Draw("E");
fLowLS->Draw("same");
auto labelLowLS = new TLatex();
labelLowLS->SetNDC();
labelLowLS->DrawLatex(
    0.53, 0.82,
    Form("#tau = %.3f #pm %.3f", fLowLS->GetParameter(1), fLowLS->GetParError(1)));
cLowLS->Draw();

std::cout << "least squares, N=50: tau = "
          << fLowLS->GetParameter(1) << " +/- " << fLowLS->GetParError(1)
          << ", status = " << static_cast<int>(resultLowLS)
          << ", ndf = " << fLowLS->GetNDF() << std::endl;

least squares, N=50: tau = 0.7887 +/- 0.1905, status = 0, ndf = 7

High-statistics comparison

With 10,000 events, most bins that determine the fit contain many counts. The Poisson distribution is then close to a Gaussian, so binned Poisson likelihood and least squares give nearly the same lifetime. Agreement in this limit explains why least squares is often effective; it does not justify using it automatically for sparse count data.

h_high_poisson = h_high.Clone("h_high_poisson_py")
h_high_ls = h_high.Clone("h_high_ls_py")
f_high_poisson = ROOT.TF1("f_high_poisson_py", "[0]*exp(-x/[1])", 0, 10)
f_high_ls = ROOT.TF1("f_high_ls_py", "[0]*exp(-x/[1])", 0, 10)
for function in (f_high_poisson, f_high_ls):
    function.SetParNames("A", "tau")
    function.SetParameters(h_high.GetMaximum(), 1.0)
    function.SetParLimits(0, 0, 1.0e9)
    function.SetParLimits(1, 0.05, 5.0)
f_high_poisson.SetLineColor(ROOT.kRed)
f_high_ls.SetLineColor(ROOT.kBlue + 1)
f_high_ls.SetLineStyle(2)

result_high_poisson = h_high_poisson.Fit(f_high_poisson, "LIRSQ0")
result_high_ls = h_high_ls.Fit(f_high_ls, "IRSQ0")

c_high = ROOT.TCanvas("c_high_py", "High-statistics comparison", 700, 500)
c_high.SetLogy()
h_high_poisson.SetMinimum(0.3)
h_high_poisson.SetTitle("10,000 events;decay time t;counts / bin")
h_high_poisson.Draw("E")
f_high_poisson.Draw("same")
f_high_ls.Draw("same")
legend_high = ROOT.TLegend(0.51, 0.72, 0.88, 0.88)
legend_high.AddEntry(
    f_high_poisson,
    f"Poisson: #tau = {f_high_poisson.GetParameter(1):.3f} #pm {f_high_poisson.GetParError(1):.3f}",
    "l",
)
legend_high.AddEntry(
    f_high_ls,
    f"least squares: #tau = {f_high_ls.GetParameter(1):.3f} #pm {f_high_ls.GetParError(1):.3f}",
    "l",
)
legend_high.Draw()
c_high.Draw()

print(
    f"binned Poisson, N=10000: tau = {f_high_poisson.GetParameter(1):.4f} "
    f"+/- {f_high_poisson.GetParError(1):.4f}, status = {int(result_high_poisson)}"
)
print(
    f"least squares, N=10000: tau = {f_high_ls.GetParameter(1):.4f} "
    f"+/- {f_high_ls.GetParError(1):.4f}, status = {int(result_high_ls)}"
)

auto hHighPoisson = static_cast<TH1D*>(hHigh->Clone("hHighPoisson"));
auto hHighLS = static_cast<TH1D*>(hHigh->Clone("hHighLS"));
auto fHighPoisson = new TF1("fHighPoisson", "[0]*exp(-x/[1])", 0, 10);
auto fHighLS = new TF1("fHighLS", "[0]*exp(-x/[1])", 0, 10);
for (auto function : {fHighPoisson, fHighLS}) {
    function->SetParNames("A", "tau");
    function->SetParameters(hHigh->GetMaximum(), 1.0);
    function->SetParLimits(0, 0, 1.0e9);
    function->SetParLimits(1, 0.05, 5.0);
}
fHighPoisson->SetLineColor(kRed);
fHighLS->SetLineColor(kBlue + 1);
fHighLS->SetLineStyle(2);

TFitResultPtr resultHighPoisson = hHighPoisson->Fit(fHighPoisson, "LIRSQ0");
TFitResultPtr resultHighLS = hHighLS->Fit(fHighLS, "IRSQ0");

auto cHigh = new TCanvas("cHigh", "High-statistics comparison", 700, 500);
cHigh->SetLogy();
hHighPoisson->SetMinimum(0.3);
hHighPoisson->SetTitle("10,000 events;decay time t;counts / bin");
hHighPoisson->Draw("E");
fHighPoisson->Draw("same");
fHighLS->Draw("same");
auto legendHigh = new TLegend(0.51, 0.72, 0.88, 0.88);
legendHigh->AddEntry(
    fHighPoisson,
    Form("Poisson: #tau = %.3f #pm %.3f", fHighPoisson->GetParameter(1), fHighPoisson->GetParError(1)),
    "l");
legendHigh->AddEntry(
    fHighLS,
    Form("least squares: #tau = %.3f #pm %.3f", fHighLS->GetParameter(1), fHighLS->GetParError(1)),
    "l");
legendHigh->Draw();
cHigh->Draw();

std::cout << "binned Poisson, N=10000: tau = "
          << fHighPoisson->GetParameter(1) << " +/- " << fHighPoisson->GetParError(1)
          << ", status = " << static_cast<int>(resultHighPoisson) << "\n"
          << "least squares, N=10000: tau = "
          << fHighLS->GetParameter(1) << " +/- " << fHighLS->GetParError(1)
          << ", status = " << static_cast<int>(resultHighLS) << std::endl;

binned Poisson, N=10000: tau = 0.9969 +/- 0.0101, status = 0
least squares, N=10000: tau = 0.9916 +/- 0.0099, status = 0

Reading the result

• The unbinned result uses every measured decay time and has no histogram-bin choice.
• The binned Poisson fit is the natural ROOT fit for a histogram of unweighted event counts, especially at low statistics.
• Least squares becomes a reasonable approximation when the fitted values and their uncertainties are approximately Gaussian; the high-statistics sample illustrates this limit.
• Use an extended likelihood when the expected total yield is a fitted part of the physics model. The RooFit examples on the reference page show that case explicitly.

In every method, inspect the fit status, parameter uncertainties, fit range, and whether the statistical model matches how the data were recorded. A small reported error is not useful if the likelihood or error model is inappropriate.