加权拟合与拟合结果的误差传播¶
以下以直线为例；由参数协方差传播拟合函数不确定度的方法同样适用于其他拟合函数。
下面的代码可在 Python / PyROOT 和 ROOT C++ 之间切换；拟合结果和图共用。
代码语言：
A weighted fit means each data point contributes to the fit according to its measurement uncertainty.
Points with smaller errors are given more weight, and points with larger errors are given less influence.
- To make sure that more precise measurements dominate the fit.
- Common in physics, engineering, and data analysis whenever error bars are known.
Mathematically, it minimizes
$$
\chi^2 = \sum_i \frac{(y_i - f(x_i))^2}{\sigma_i^2}
$$
where $\sigma_i$ is the error of each data point.
For a TGraphErrors, ROOT reads the Y errors and uses them in the fit:
auto r = graph->Fit(fitfunc, "QS");
Q suppresses ROOT's standard fit printout, while S returns the complete fit result. Do not add W here: for TGraphErrors, W ignores the point errors.
After fitting, the uncertainty of the fitted function $f(x)$ can be evaluated at any x-value:
TGraphErrors gpt(1);
gpt.SetPoint(0, x0, 0);
TVirtualFitter::GetFitter()->GetConfidenceIntervals(&gpt, 0.68);
double fx = gpt.GetY()[0]; // fitted value
double efx = gpt.GetEY()[0]; // uncertainty of f(x)
ROOT uses the covariance matrix of the fitted parameters to calculate $$ \sigma_{f(x)}^2 = J\,\mathrm{Cov}(p)\,J^T, $$ where $J$ is the gradient of $f(x)$ with respect to its parameters.
import ROOT
from array import array
x = array("d", [0, 2, 4, 6, 8, 10])
y = array("d", [5, 8, 14, 20, 24, 35])
ex = array("d", [0] * len(x))
ey = array("d", [2, 1.5, 2.5, 1, 3, 2])
# TGraphErrors supplies ey to the weighted fit. W must not be used here.
graph = ROOT.TGraphErrors(len(x), x, y, ex, ey)
graph.SetTitle("Weighted linear fit; x; y")
graph.SetMarkerStyle(20)
graph.SetMarkerColor(ROOT.kBlue + 1)
fitfunc = ROOT.TF1("fitfunc_py", "pol1", 0, 10)
fitfunc.SetParNames("c", "m")
result = graph.Fit(fitfunc, "QS")
print(f"c = {fitfunc.GetParameter(0):.4f} +/- {fitfunc.GetParError(0):.4f}")
print(f"m = {fitfunc.GetParameter(1):.4f} +/- {fitfunc.GetParError(1):.4f}")
print(f"chi2/ndf = {result.Chi2():.4f} / {result.Ndf()} = {result.Chi2()/result.Ndf():.4f}")
# GetConfidenceIntervals propagates the parameter covariance to f(x).
band = ROOT.TGraphErrors(200)
for i in range(200):
band.SetPoint(i, 10.0 * i / 199.0, 0.0)
ROOT.TVirtualFitter.GetFitter().GetConfidenceIntervals(band, 0.68)
x0 = 5.0
point = ROOT.TGraphErrors(1)
point.SetPoint(0, x0, 0.0)
ROOT.TVirtualFitter.GetFitter().GetConfidenceIntervals(point, 0.68)
fx0, efx0 = point.GetY()[0], point.GetEY()[0]
print(f"f({x0:.1f}) = {fx0:.4f} +/- {efx0:.4f} (68% confidence interval)")
canvas = ROOT.TCanvas("c_weighted_py", "Weighted fit", 800, 600)
graph.Draw("AP")
band.SetFillColor(ROOT.kOrange)
band.SetFillStyle(3001)
band.Draw("3 same")
fitfunc.SetLineColor(ROOT.kGreen + 2)
fitfunc.SetLineWidth(2)
fitfunc.Draw("same")
point.SetMarkerStyle(21)
point.SetMarkerColor(ROOT.kRed + 1)
point.SetLineColor(ROOT.kRed + 1)
point.Draw("PE same")
legend = ROOT.TLegend(0.12, 0.70, 0.42, 0.90)
legend.AddEntry(graph, "Data", "PE")
legend.AddEntry(fitfunc, "Fit (pol1)", "L")
legend.AddEntry(band, "68% confidence band", "F")
legend.AddEntry(point, "f(5) with propagated error", "PE")
legend.Draw()
box = ROOT.TPaveText(0.53, 0.68, 0.89, 0.90, "NDC")
box.SetFillColor(ROOT.kWhite)
box.SetBorderSize(1)
box.SetTextAlign(12)
box.AddText(f"c = {fitfunc.GetParameter(0):.3f} #pm {fitfunc.GetParError(0):.3f}")
box.AddText(f"m = {fitfunc.GetParameter(1):.3f} #pm {fitfunc.GetParError(1):.3f}")
box.AddText(f"#chi^{{2}}/ndf = {result.Chi2():.3f}/{result.Ndf()} = {result.Chi2()/result.Ndf():.3f}")
box.AddText(f"f(5) = {fx0:.3f} #pm {efx0:.3f}")
box.Draw()
canvas.Draw()
In [1]:
#include "TGraphErrors.h"
#include "TF1.h"
#include "TVirtualFitter.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TPaveText.h"
#include "TString.h"
#include <iostream>
#include <iomanip>
// Data and their Y uncertainties
double x[] = {0, 2, 4, 6, 8, 10};
double y[] = {5, 8, 14, 20, 24, 35};
double ey[] = {2, 1.5, 2.5, 1, 3, 2};
const int n = sizeof(x) / sizeof(double);
// TGraphErrors stores the measured values and their uncertainties.
auto graph = new TGraphErrors(n, x, y, nullptr, ey);
graph->SetTitle("Weighted linear fit; x; y");
graph->SetMarkerStyle(20);
graph->SetMarkerColor(kBlue + 1);
double xmin = 0.0, xmax = 10.0;
auto fitfunc = new TF1("fitfunc", "pol1", xmin, xmax);
fitfunc->SetParNames("c", "m");
// Q: quiet output; S: return TFitResult. TGraphErrors supplies the weights through ey.
auto r = graph->Fit(fitfunc, "QS");
std::cout << std::fixed << std::setprecision(4);
std::cout << "Fit result\n";
std::cout << "c = " << fitfunc->GetParameter(0)
<< " +/- " << fitfunc->GetParError(0) << "\n";
std::cout << "m = " << fitfunc->GetParameter(1)
<< " +/- " << fitfunc->GetParError(1) << "\n";
std::cout << "chi2/ndf = " << r->Chi2() << " / " << r->Ndf()
<< " = " << r->Chi2() / r->Ndf() << std::endl;
Fit result c = 2.7285 +/- 1.3242 m = 2.9692 +/- 0.2316 chi2/ndf = 4.1898 / 4 = 1.0475
In [2]:
// Evaluate the 68% confidence band of the fitted mean.
const int nb = 200;
auto gMean = new TGraphErrors(nb);
for (int i = 0; i < nb; ++i) {
double xi = xmin + (xmax - xmin) * i / (nb - 1.0);
gMean->SetPoint(i, xi, 0.0);
}
TVirtualFitter::GetFitter()->GetConfidenceIntervals(gMean, 0.68);
// Evaluate the fitted value and its propagated uncertainty at x = 5.
double x0 = 5.0;
TGraphErrors gx0(1);
gx0.SetPoint(0, x0, 0.0);
TVirtualFitter::GetFitter()->GetConfidenceIntervals(&gx0, 0.68);
double fx0 = gx0.GetY()[0];
double efx0 = gx0.GetEY()[0];
std::cout << std::fixed << std::setprecision(4);
std::cout << "f(" << x0 << ") = " << fx0
<< " +/- " << efx0 << " (68% confidence interval)" << std::endl;
f(5.0000) = 17.5745 +/- 0.6683 (68% confidence interval)
In [3]:
auto c1 = new TCanvas("c1", "Weighted fit", 800, 600);
graph->Draw("AP");
// Draw the confidence band behind the fitted function.
gMean->SetFillColor(kOrange);
gMean->SetFillStyle(3001);
gMean->Draw("3 same");
fitfunc->SetLineColor(kGreen + 2);
fitfunc->SetLineWidth(2);
fitfunc->Draw("same");
// Mark f(5) and its propagated uncertainty.
auto px0 = new TGraphErrors(1);
px0->SetPoint(0, x0, fx0);
px0->SetPointError(0, 0.0, efx0);
px0->SetMarkerStyle(21);
px0->SetMarkerColor(kRed + 1);
px0->SetLineColor(kRed + 1);
px0->Draw("PE same");
auto leg1 = new TLegend(0.12, 0.70, 0.42, 0.90);
leg1->AddEntry(graph, "Data", "PE");
leg1->AddEntry(fitfunc, "Fit (pol1)", "L");
leg1->AddEntry(gMean, "68% confidence band", "F");
leg1->AddEntry(px0, "f(5) with propagated error", "PE");
leg1->Draw();
// Print the numerical fit result directly on the plot.
auto resultBox = new TPaveText(0.53, 0.68, 0.89, 0.90, "NDC");
resultBox->SetFillColor(kWhite);
resultBox->SetBorderSize(1);
resultBox->SetTextAlign(12);
resultBox->AddText(Form("c = %.3f #pm %.3f", fitfunc->GetParameter(0), fitfunc->GetParError(0)));
resultBox->AddText(Form("m = %.3f #pm %.3f", fitfunc->GetParameter(1), fitfunc->GetParError(1)));
resultBox->AddText(Form("#chi^{2}/ndf = %.3f/%d = %.3f", r->Chi2(), r->Ndf(), r->Chi2() / r->Ndf()));
resultBox->AddText(Form("f(5) = %.3f #pm %.3f", fx0, efx0));
resultBox->Draw();
c1->Update();