Re: extraction of likelhood value

Hello Jan,

  you ned to replace the line:

> TVirtualFitter *fitter = TVirtualFitter::Fitter(hWimp);

by

TVirtualFitter *fitter = TVirtualFitter::GetFitter();

In the first case you were creating a new instance of TVirtualFitter, while you need to get the instance used in the last TH1::Fit to retrieve the information like the minimum of the likelihood.

Cheers

  Lorenzo

On Jul 12, 2007, at 10:34 AM, Jan Conrad wrote:

> Dear Rooters,
>
>
> I am running the included macro. The goal is to extract the
> likelihood
> value (or log likelihood, doesn't matter).
>
> If fitting using chi2, using:
>
> TVirtualFitter *fitter = TVirtualFitter::Fitter(hWimp);
> Double_t amin,edm,errdef;
> Int_t nvpar,nparx;
> fitter->GetStats(amin,edm,errdef,nvpar,nparx);
> Double_t chi2 = hWimp->GetFunction("fitSignal")->GetChisquare();
> printf("amin=%g, chi2=%g\n", amin, chi2);
>
> I get:
> amin=51.689, chi2=51.689
>
> .. which makes sense.
>
>
> If fitting using LL (log likelihood), I get:
>
>> amin=502.46, chi2=857560
>
> The obtained fit parameters are similar (though the likelihood fit
> did not
> catastrophical fail), I would have expected the chi2 to be similar
> to the
> one obtained in the first case. Furthermore in
>
> http://root.cern.ch/root/roottalk/roottalk00/1963.html
>
> Rene says that amin = -2 logL. Is that really true ?
>
> I am using ROOT v5.10/00 on
> Scientific Linux CERN SLC release 4.2 (Beryllium)
> Linux 2.6.9-22.0.1.EL.cern
>
> cheers, thanks,
> Jan
>
>
>
> ToyLL(){
> //////////////////////////////////////////////////////////
> //
> // Jan Conrad (Stockholm University), conrad_at_physto.se
> //
> //////////////////////////////////////////////////////////
>
>
> Int_t nBinsChi2 = 50;
>
>
> // generate histogram with E-2 spectrum
>
>
> TH1F* hWimp = new TH1F("wimp","wimp",nBinsChi2,1e4,1e6);
>
>
>
> // Powerlaw (function 1)
> TF1 *f1 = new TF1("f1","(x/1.e5)**(-2.0)",1e4.,1e6);
>
> // Double exponential (function 2).
> TF1 *f2 = new TF1("f2","TMath::Exp(-1.e-5*x)",1e4,1e6);
>
>
>
>
> Double_t WimpSpec(Double_t *x, Double_t *par) {
> return par[0]*TMath::Exp(par[1]*x[0]);
> }
>
>
>
>
> hWimp->Reset();
>
>
>
> for (Int_t iExp=0; iExp<1000000; iExp++) {
> double ws = f2->GetRandom();
> hWimp->Fill(ws);
> }
>
> hWimp->Sumw2();
>
>
> // check signal fitting
> TF1 *fitSignal = new TF1("fitSignal",WimpSpec,1e4.,1e6,2);
> fitSignal->SetNpx(500);
> fitSignal->SetLineWidth(4);
> fitSignal->SetLineColor(kMagenta);
> fitSignal->SetParameters(1.e4,-1.e-4);
> //TF1 *f2 = new
> TF1("f2","TMath::Exp(-1.e-4*x)+.2*TMath::Exp(-1.e-5*x)",1e4,1e6);
> //return par[0]*TMath::Exp(par[1]*x[0])+par[2]*TMath::Exp(par[3]*x
> [0]);
> hWimp->Fit("fitSignal","L,V,I");
>
>
> TVirtualFitter *fitter = TVirtualFitter::Fitter(hWimp);
> Double_t amin,edm,errdef;
> Int_t nvpar,nparx;
> fitter->GetStats(amin,edm,errdef,nvpar,nparx);
> Double_t chi2 = hWimp->GetFunction("fitSignal")->GetChisquare();
> printf("amin=%g, chi2=%g\n", amin, chi2);
> }
> ______________________________________________________________________
> __
> J. Conrad, Physics Department, Stockholm University,
> AlbaNova, SE-10691 Stockholm, Sweden, http://www.particle.kth.se/
> ~conrad
>
Received on Wed Jul 18 2007 - 17:25:17 CEST