Re: Problems fitting

From: Rene Brun (Rene.Brun@cern.ch)
Date: Fri Mar 06 1998 - 17:58:04 MET


Chris,
In case the option "R" in TH1::Fit is specified, I have discovered,
thanks to your example, that this function uses the overflow bin
in the fit!
I have fixed this problem in my development version.
Meanwhile, remove the option "R" and you will get a correct fit
and correct chisquare.

Rene Brun


Chris Jillings wrote:
> 
> Hi Rooters,
>     I have found a pathological fit and I am trying to understand why root
> failed to fit properly.
>     I have a macro below that starts the fit at a particular set of
> parameter values close to the correct answer. I calcuated the chi square
> two ways. First I looped over the histogram bins and calculated it myself
> in the macro (Answer: 5233). Second, I set very tight bounds on the
> paremeters, effectively forcing migrad not to change them, and then looked
> at the value of FCN from migrad( Answer 5E8!).
> 
> The output of the macro is
> (bin #)  data           model           (model-data)**2  total chiSq
>                                         ---------------
>                                            data
> 
>                   :
>                   :
> 
> 70      2138.000000     2790.239235     198.978494      1390.020928
> 71      2547.000000     3419.200307     298.678200      1688.699128
> 72      3267.000000     4196.629537     264.527418      1953.226546
> 73      4010.000000     5158.748900     329.083300      2282.309846
> 74      4942.000000     6350.814741     401.610477      2683.920323
> 75      6087.000000     7829.403125     498.762715      3182.683039
> 76      7695.000000     9665.278362     504.483018      3687.166056
> 77      9704.000000     11946.994351    518.448440      4205.614496
> 78      12155.000000    14785.416850    569.238404      4774.852900
> 79      15641.000000    18319.403236    458.656345      5233.509246
>  **********
>  **    1 **SET ERR           1
>  **********
>  FCN=5.16288e+08 FROM MIGRAD    STATUS=CONVERGED      48 CALLS          49
> TOTAL
>                      EDM=1.25522e-06    STRATEGY= 1      ERROR MATRIX
> ACCURATE
>   EXT PARAMETER                                   STEP         FIRST
>   NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE
>    1  1st slope    9.20000e+00   9.23121e-08   2.13212e-02** at limit **
>    2  2nd Amp      2.00000e-01   1.20541e-09   7.70472e-03** at limit **
>    3  2nd slope    5.70000e+00   4.82292e-07   4.87308e-02** at limit **
>    4  scale fac    1.70000e+04   1.64631e-05   9.00421e-03** at limit **
> 
> 
> So FCN, which from the docs is a chi square, is 5E8, whereas my
> loop calculates 5233.
> 
> Any ideas? Thanks once again for your help. The full macro follows.
> 
> Chris
> 
>  {
>   gROOT->Reset();
>   gStyle->SetOptFit(1);
>   gStyle->SetStatX(0.45);
>   gStyle->SetStatY(0.85);
> 
>   TCanvas* c1 = new TCanvas("c1","Ntuple Plots",1); // create canvas
>   c1->SetFillColor(kWhite);
> 
>   TFile* f1 = new TFile("cjj_g10.root","READ");
>   TH1F* hi1 = new TH1F("hi1","Angular Resolution: Gamma Energy = 10 MeV
> ",80,-1.0,1.0);
>   hi1->SetXTitle("Cos`q#");
>   //  h509->Draw("Angres>>hi1","Egen>4&&Rfitt<600","goff");
>   h503->Draw("(Ue*Uf+Ve*Vf+We*Wf)>>hi1","Rfit<600","goff");
>   c1->SetLogy(1);
>   TF1* expExp = new TF1("expExp",expExp,-1,1,4);
>   expExp->SetParameters(9.2,0.2,5.7,17000);
>   expExp->SetParNames("1st slope","2nd Amp","2nd slope","scale fac");
>   Double_t x;
>   Double_t y;
>   Double_t dChiSq;
>   Double_t chiSq = 0..0;
>   Double_t data;
>   Double_t sigSq;
>   printf("Hello\n");
>   for( Int_t i=0 ; i<80 ; i++) {
>     x = -0.9875+i*0.025;
>     //    printf("%d\t%f\n",i,x);
>     y = expExp->Eval(x);
>     data = hi1->GetBinContent(i);
>     if (data==0.0) sigSq = y;
>     else sigSq = data;
>     dChiSq = (data-y)*(data-y)/sigSq;
>     chiSq += dChiSq;
>     printf("%d\t%f\t%f\t%f\t%f\n",i,data,y,dChiSq,chiSq);
>   }
>   //  expExp->Draw();
>   expExp->SetParLimits(0,9.2,9.3);
>   expExp->SetParLimits(1,0.2,0.21);
>   expExp->SetParLimits(2,5.7,5.8);
>   expExp->SetParLimits(3,17000.0,17100.0);
> 
>   TPostScript ps("cjj_g10.eps",114);
>   hi1->Fit("expExp","RB");
>   //  expExp->Draw();
>   //  hi1->Draw("same");
>   c1->Update();
>   ps.Close();
> 
> }



This archive was generated by hypermail 2b29 : Tue Jan 04 2000 - 00:34:30 MET