Hi ROOTers,
I would like to fit the data using lognormal distributiuon and
I wrote attached macro.
But, I alwasy got failure until I try to fit 3 times.
As shown below, in 3 times, I got something good results.
However, if I change value of npars, which is argument of macro,
for example, to 4 or 2, I succeeded in fitting the histogram a time.
Why ?
I use ROOT v 3.03/05 and g++ 2.95.2 on IRIX 64.
Harufumi Tsuchiya
--------------------------------------------------------
root [0] .x LogNormal2.cpp() <------ (first try)
200
FCN=8096.26 FROM MIGRAD STATUS=CONVERGED 555 CALLS 556 TOTAL
EDM=2.13927e-07 STRATEGY= 1 ERR MATRIX NOT POS-DEF
EXT PARAMETER APPROXIMATE STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 p0 1.02155e+01 3.06869e+00 1.02872e-02 -2.64182e-03
2 p1 8.28982e+05 3.63901e+07 2.54515e+02 2.29469e-10
3 p2 9.99117e+06 9.57300e+07 2.54515e+02 -1.94306e-11
Chi2 and dof = 8096.26 184 <---- wrong !!!!!!!!!!!!!!!!!!!!!!
Probability 0
root [1] .x LogNormal2.cpp() <------- (second try)
Warning in <TH1::Build>: Replacing existing histogram: hl (Potential memory leak).
Warning in <TH1::Build>: Replacing existing histogram: hl2 (Potential memory leak).
200
FCN=8019.71 FROM MIGRAD STATUS=CONVERGED 508 CALLS 509 TOTAL
EDM=1.11024e-07 STRATEGY= 1 ERR MATRIX NOT POS-DEF
EXT PARAMETER APPROXIMATE STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 p0 1.10649e+01 2.78345e+00 1.08733e-02 -1.87808e-03
2 p1 8.72864e+06 4.67588e+07 2.53309e+02 1.32604e-10
3 p2 3.70683e+07 2.19576e+08 2.53309e+02 -3.13716e-11
Chi2 and dof = 8019.71 181 <---- wrong !!!!!!!!!!!!!!!!!!!!!!
Probability 0
root [2] .x LogNormal2.cpp() <----- (third try)
Warning in <TH1::Build>: Replacing existing histogram: hl (Potential memory leak).
Warning in <TH1::Build>: Replacing existing histogram: hl2 (Potential memory leak).
200
FCN=199.661 FROM MIGRAD STATUS=CONVERGED 469 CALLS 470 TOTAL
EDM=1.99527e-07 STRATEGY= 1 ERROR MATRIX UNCERTAINTY 2.2 per cent
EXT PARAMETER STEP FIRST
NO. NAME VALUE ERROR SIZE DERIVATIVE
1 p0 1.18611e+02 1.45605e+00 1.39037e-02 -1.25445e-04
2 p1 1.49387e+00 5.31796e-03 -7.32891e-06 -1.45410e-01
3 p2 4.96796e-01 3.59119e-03 -3.15101e-05 -9.75802e-02
Chi2 and dof = 199.661 185 <---- something good !!
Probability 0.21847
-----------------------------------------------------------------------------------
void LogNormal2(const Int_t npars = 3) {
const Int_t nbins = 200;
Double_t x[npars], q;
Double_t lx[npars], lq;
Float_t xmin, xmax;
xmin = 0; xmax = npars;
TH1F *hl = new TH1F("hl", "hl", nbins, xmin, xmax);
Int_t ndata = 0;
Float_t xbins[nbins+1];
for ( Int_t i = 0; i < nbins; i++ ) {
xbins[ndata] = xmin + (xmax-xmin)/nbins*Float_t(i);
xbins[ndata] = 10**xbins[ndata];
ndata++;
}
TH1F *hl2 = new TH1F("hl2", "hl2", ndata-1, xbins);
cout << ndata << endl;
lq = 0.0;
for ( Int_t i = 0; i < 10000; i++) {
for ( Int_t j = 0; j < npars; j++ ) {
lx[j] = gRandom->Rndm(i); x[j] = 10**lx[j];
lq += lx[j];
}
q = 10**lq;
// cout << q << " " << lq << endl;
hl->Fill(lq); hl2->Fill(q);
lq = 0.0;
}
gROOT->SetStyle("Plain");
gStyle->SetOptTitle(0);
TCanvas *c1 = new TCanvas("c1", "", 10, 10, 500, 500);
hl->Draw();
TCanvas *c2 = new TCanvas("c2", "", 300, 10, 500, 500);
c2->SetLogx();
hl2->Draw();
// Define log gauss //
Float_t l, h;
l = 10**xmin; h = 10**xmax;
TF1 *f1 = new TF1("f1", "[0]*exp(-(log10(x) - [1])**2/(2*[2]**2))", l, h);
Float_t par[3];
par[0] = 200;
par[1] = TMath::Log10(300.0);
par[2] = TMath::Log(100);
f1->SetParameters(par[0], par[1], par[2]);
hl2->Fit("f1", "R");
Double_t chi2;
Int_t dof;
chi2 = f1->GetChisquare(); dof = f1->GetNDF();
cout << "Chi2 and dof = " << chi2 << " " << dof << endl;
cout << "Probability " << TMath::Prob(chi2, dof) << endl;
}
This archive was generated by hypermail 2b29 : Sat Jan 04 2003 - 23:50:57 MET