[ROOT] Problems with fitting with chi-square

From: Alexander Dietz (Alexander.Dietz@mpi-hd.mpg.de)
Date: Tue Mar 26 2002 - 11:30:09 MET 

here,

I have a problem when fitting with the chi2-square method!

With the following program I create several random-histograms consisting
of two parts: a flat poisson-background (with mean 2) and a Gaussian line
centered at channel 50, with sigma 2.0 and with 100 counts in it.
Then I fit this histogram with the function 'function' (assuming a flat
backgriund and a gaussian line), using the Chi-square method or the
Maximum likelihood method (adding a 'L' to the fitting options).
Each results of each fit are then given into a second histogram 'sum'.

When I do the fit with Maximum likelihood, the results in histogram 'sum'
are really good, fluctuating around the true value of 100. But when I use
the chi-square method, the results are centered around 95 which is quite
wrong! Only a few percent of the results give values larger than 100!

So what is wrong here? Maybe the chi-square method itself gives confusing
results? Or is there a bug in the fitting procedure?

Thanks for help

	A. Dietz


P.S. When doing the same but with 10000 counts in the line, BOTH methods
give very good results.


Program:


start()
{
  gROOT->Reset();

  c1 = new TCanvas("c1","The fitting problem",200,10,600,400);
  c1->SetGrid();


  // Create one histograms.
  main   = new TH1F("main","Test spectrum",100,0,100);
  sum    = new TH1F("sum" ,"Fittig results",100,80,120);
  sum->SetStats(kFALSE);
  sum->SetFillColor(47);
  main->SetStats(kFALSE);
  main->SetFillColor(46);

  for (Int_t k=0;k<=100;k++) {

    // setting random seeds
    gRandom->SetSeed(k);

    Float_t value;
    // filling spectrum with Poisson background
    for (Int_t i=0;i<=100;i++) {
      value=gRandom->Poisson(1);
      main->SetBinContent(i,value);
    }

    // adding gaussian line at channel 50, sigma 2.0=channel with 100
counts in it
    Int_t counts=100;
    for ( Int_t i=0; i<counts; i++) {
      value = gRandom->Gaus(50,2.0);
      main->Fill(value);
    }

    main->Draw("");
    c1->Update();
    c1->Modified();

    // fiting the spectrum
    TF1* func=new TF1("func",function,0,100,4);
    func->SetParameters(1,100,50,2);
    main->Fit("func","QWEMR"); // or with likelihood: "LQWEMR"

    // printing out size of line
    Float_t size=func->GetParameter(1);
    cout << "Size of line: " << size << endl;

    // store result to histogram 'sum'
    sum->Fill(size);
  }

  // showing the results of the fit
  sum->Draw();
}

// fitting function
Double_t function(Double_t* x, Double_t* p)
{
  Double_t dummy=(x[0]-p[2])/p[3];
  Double_t result=p[0] + p[1]*exp(-0.5*dummy*dummy)/ (2.506 * p[3]);
  return result;
}

This archive was generated by hypermail 2b29 : Sat Jan 04 2003 - 23:50:47 MET