Fitten korrelierter Daten

Next: Güte des Fit - Up: Fitten mit ROOT Previous: Unbinned Maximum-Likelihood Fit

Fitten korrelierter Daten

Die einfache Formel

$\ensuremath{\displaystyle \ensuremath{\chi^2}\equiv \sum_{i=1}^N \frac{( y_i - f(x_i,\vec{a} ))^2} {(\sigma_i)^2}}$

gilt nur für unkorellierte Messwerte, es gibt aber häufig Korrelationen zwischen Messwerten, z.B.

gemeinsame Normierungsunsicherheit
gemeinsamer Nullpunkt/Startwert
andere systematische Effekte

Dann muss die `Kovarianzmatrix' für den Fit benutzt werden. Für $\ensuremath{\displaystyle n}$ Werte $\ensuremath{\displaystyle y_1, y_2, ..., y_n}$ beschreibt die $\ensuremath{\displaystyle n\times n}$ Kovarianzmatrix $\ensuremath{\displaystyle \cal{C}}$ Fehler und Korrelationen der einzelnen Werte:

$\begin{displaymath} {\cal{C}} = \left( \begin{array}{cccc} \sigma_{y_1}^2 & cov... ...& cov_{y_2,y_n} & ... & \sigma_{y_n}^2 \\ \end{array} \right) \end{displaymath}$

mit $\ensuremath{\displaystyle cov_{y_i,y_j} = \rho_{ij} \sigma_{y_i} \sigma_{y_j}}$ .

Das $\ensuremath{\displaystyle \ensuremath{\chi^2}}$ wird dann gemäss der allgemeinen Formel

$\ensuremath{\displaystyle \ensuremath{\chi^2}\equiv \vec{\Delta}^{ T} \cdot {\cal{C}}^{-1} \cdot \vec{\Delta}}$

berechnet werden
( $\ensuremath{\displaystyle {\cal{C}}^{-1} = \mbox{invertierte Kovarianzmatrix}}$ )
( $\ensuremath{\displaystyle \vec{\Delta} = \mbox{ Vektor der Residuen :} ( y_i - f(x_i, \vec{a}) )}$ )

Zum Erstellen und Invertieren der Kovarianzmatrix $\ensuremath{\displaystyle {\cal{C}}}$ empfiehlt sich die ROOT Klasse TMatrixD.

Resonanzkurve des Wirkungsquerschnitts

Gutes Beispiel für korrelierte Fehler ist Messung des Wirkungsquerschnitts bei LEP im Bereich der Resonanz.

$\textstyle \parbox{.5\textwidth}{ Experimentell: \begin{center} {\ensuremath{\... ...in $x$ aus Unsicherheit in {\ensuremath{\displaystyle E_{CM}}}\end{itemize}}$ $\textstyle \parbox{.45\textwidth}{ \begin{center} \mbox{ \includegraphics[width=.45\textwidth]{ls.eps}} \end{center} }$

Klassische Messung, auch im Fortgeschrittenen Praktikum an der LMU !

$\textstyle \parbox{1.5cm}{ {\color{blue}$e^+e^-$} (3.3 \%) }$ $\textstyle \parbox{10.1cm}{ \includegraphics[height=5cm,bb=19 24 566 572]{ee.eps} \includegraphics[height=5cm,bb=19 24 566 572]{mm.eps} }$ $\textstyle \parbox{1.2cm}{ {\color{blue} $\mu^+\mu^-$} (3.3 \%)}$
$\textstyle \parbox{1.5cm}{ {\color{blue}$\tau^+\tau^-$} (3.3 \%) }$ $\textstyle \parbox{10.1cm}{ \includegraphics[height=5cm,bb=19 24 566 572]{tt.eps} \includegraphics[height=5cm,bb=19 24 566 572]{mh.eps} }$ $\textstyle \parbox{1.2cm}{ {\color{blue} $q\bar{q}$} (70 \%)}$

Entsprechender Fit mit Minuit:

#include "TMatrix.h"

TMatrixD *cov;
double *val;
double *xval;
const int nv = 7;

void setInput() 
{ // define input
  val = new double[nv];
  xval = new double[nv];
  cov = new TMatrixD(nv,nv);

  double xv[nv] = { 88.396, 89.376, 90.234, 91.238, 92.068, 93.080, 93.912 };
  double yv[nv] = {  6.943, 13.183, 25.724, 40.724, 27.031, 12.273, 6.980 };
  double ey[nv] = {  0.087, 0.119, 0.178, 0.087, 0.159, 0.095, 0.064};
  double rsys = 0.004; // correllated relative syst uncertainty on yv

  for ( int i=0; i<nv; i++ ) {
    val[i] = yv[i];
    xval[i] = xv[i];
    cov(i,i) = ey[i]*ey[i]; // stat errors
    for ( int j=i; j<nv; j++ ) {
      cov(i,j) = cov(i,j) + yv[i] * yv[j] * rsys * rsys;
      cov(j,i) = cov(i,j);
    }
  }
  
  cov->Print();
  cov->Invert();
  cov->Print();
}

double BreitWig( double *x, double *par) 
{
  // Breit- Wigner function
  double mw = par[0], gw = par[1], spk=par[2], mw2, gw2, eb2;
  
  mw2 = mw*mw;
  gw2 = gw*gw;
  eb2 = x[0]*x[0];
  return( spk * gw2*mw2 / ( pow( eb2 - mw2, 2 ) + mw2 * gw2 ) );
  //  return( spk * gw2*mw2 / ( pow( eb2 - mw2, 2 ) + eb2*eb2/mw2 * gw2 ) );
}

void fcn(Int_t &npar,double *gin, double &f, double *par, int iflag)
 {
   // calculate Chi2
   double chi2 = 0.;
   for ( int i=0; i<nv; i++ ) {
     double fvali = BreitWig( &xval[i], par );
     for ( int j=0; j<nv; j++ ) {
       double fvalj = BreitWig( &xval[j], par );
       chi2 += ( fvali - val[i] ) * cov(i,j) * ( fvalj - val[j] ) ;
     }
   }
   //   cout << "Chi2 = " << chi2 << endl;
   f = chi2;

 }

void bwfitmnc()
{

  setInput();

  //   initialize TMinuit 
  const int npar = 3;
  gMinuit = new TMinuit(npar);  
  // Tell Minuit the function 
  gMinuit->SetFCN(fcn);
  // Start Werte fuer Minuit
  int ierflg;
  gMinuit->mnparm(0, "mass", 85, 0.1, 0., 0., ierflg);
  gMinuit->mnparm(1, "width", 2., 0.1, 0., 0., ierflg);
  gMinuit->mnparm(2, "speak", 10., 0.1, 0., 0., ierflg);
  // minimization
  double arglist[3] = { 1000., 1., 0. };
  gMinuit->mnexcm("MIGRAD", arglist ,1,ierflg);
  
  //  fcn( nparfit, &zero, &chi2, par, 5, 0 );

   //  TGraphErrors *tg = new TGraphErrors( 7, xv, yv, ex, ey );
   //tg->Draw("AP");
  // central value
  double par[npar], epar[npar];
  for ( int i=0; i<npar; i++ ) {
    gMinuit->GetParameter( i, par[i], epar[i] ) ;

  }

  double ey[nv] = {  0.087, 0.119, 0.178, 0.087, 0.159, 0.095, 0.064};
  for ( int i=0; i<nv; i++ ) {
    //    double delta = (yv[i]  - f->Eval(xv[i]);
    double delta = (val[i]  - BreitWig(&xval[i], par));
    cout << val[i] << " +- " << ey[i]  << " delta: " <<  delta << "   " << delta/ey[i] << endl;
  }
}

Berücksichtigung von Fehlern in der x-Koordinate

Bei Messungen von $\ensuremath{\displaystyle (x,y)}$ Wertepaaren gibt es oft Fehler nicht nur in $\ensuremath{\displaystyle \Delta y}$ sondern auch in $\ensuremath{\displaystyle \Delta x}$ .

Der übliche $\ensuremath{\displaystyle \ensuremath{\chi^2}}$ Ansatz erlaubt nur Unsicherheiten/Kovarianz in $\ensuremath{\displaystyle y}$
Trick: Über Steigung der Fit-Funktion kann man Fehler $\ensuremath{\displaystyle \Delta x}$ in Fehler $\ensuremath{\displaystyle \Delta y}$ transformieren:
$\ensuremath{\displaystyle \sigma_i^{y(x)} = \left. \frac{d f(x, \vec{a})}{ d x}\right\vert _{x=x_i} \sigma_i^{x} }$
Iterativ vorgehen:
- Erst Fit ohne Berücksichtigung der $\ensuremath{\displaystyle \Delta x}$ $\ensuremath{\color{dgreen}{\Rightarrow}}$ Näherung für Fit-Funktion und Parameter
- Dann Steigung der Fit-Funktion an den Stellen $\ensuremath{\displaystyle x_i}$ bestimmen und Fehler zur Kovarianzmatrix addieren.
- Meist genügt einer Iteration ...

#include "TMatrix.h"

TMatrixD *cov;
double *val;
double *xval;
const int nv = 7;

void setInput( bool deltax = true, double * par = 0 ) 
{ // define input
  val = new double[nv];
  xval = new double[nv];
  cov = new TMatrixD(nv,nv);

  double xv[nv] = { 88.396, 89.376, 90.234, 91.238, 92.068, 93.080, 93.912 };
  double yv[nv] = {  6.943, 13.183, 25.724, 40.724, 27.031, 12.273, 6.980 };
  double ey[nv] = {  0.087, 0.119, 0.178, 0.087, 0.159, 0.095, 0.064};
  double rsys = 0.004; // correlated relative syst uncertainty on yv
  double esys = 0.007; // 7 GeV syst error on x=E (absolute)

  TF1 * f1;

  if ( deltax ) {
    // use as Root function
    TF1 * f1 = new TF1( "bw", BreitWig, 70, 100, 3 );
    f1->SetParameters( par );
  }
  for ( int i=0; i<nv; i++ ) {
    val[i] = yv[i];
    xval[i] = xv[i];
    (*cov)(i,i) = ey[i]*ey[i]; // stat errors
    for ( int j=i; j<nv; j++ ) {
      (*cov)(i,j) = (*cov)(i,j) + yv[i] * yv[j] * rsys * rsys;

      if ( deltax ) { // treat errors in x, use TF1 Derivative method
    double dydxi = f1->Derivative( xv[i] ); // slope at x[i]
    double dydxj = f1->Derivative( xv[j] ); // slope at x[j]
    (*cov)(i,j) = (*cov)(i,j) + dydxi * esys * dydxj * esys;
      }

      (*cov)(j,i) = (*cov)(i,j);
    }
  }
  cov->Print();
  cov->Invert();
  cov->Print();
}
...
void bwfitmnce()
{

  setInput(false);
  //   initialize TMinuit 
  const int npar = 3;
  gMinuit = new TMinuit(npar);  
  // Tell Minuit the function 
  gMinuit->SetFCN(fcn);
  // Start Werte fuer Minuit
  int ierflg;
  gMinuit->mnparm(0, "mass", 85, 0.1, 0., 0., ierflg);
  gMinuit->mnparm(1, "width", 2., 0.1, 0., 0., ierflg);
  gMinuit->mnparm(2, "speak", 10., 0.1, 0., 0., ierflg);
  // minimization
  double arglist[3] = { 1000., 1., 0. };
  gMinuit->mnexcm("MIGRAD", arglist ,1,ierflg);
  
  double par[npar], epar[npar];
  for ( int i=0; i<npar; i++ ) {
    gMinuit->GetParameter( i, par[i], epar[i] ) ;
  }

  // include x errors and iterate fit 
  setInput(true, par);
  gMinuit->mnexcm("MIGRAD", arglist ,1,ierflg);
  // accurate error matrix
  gMinuit->mnexcm("HESSE", arglist ,1,ierflg);

  ...  
}

Next: Güte des Fit - Up: Fitten mit ROOT Previous: Unbinned Maximum-Likelihood Fit

GDuckeck 2008-07-01