Hi Hajime,
Once more you are correct. Fix now in CVS.
Rene Brun
On Sat, 11 Sep 2004, Hajime
Nanjyo wrote:
> Dear Rene,
>
> I think "logL += TMath::Log(prob) * weight[i]/sum;" should be
> "logL += TMath::Log(prob) * weight[i];".
>
> Don't divide "weight[i]" by "sum". I'm sorry .
> This makes difference in the error evaluation of the fit.
> I think the fix below is right...
>
> //////////////////////////////////////////////////////////
> prob = fitfunc->EvalPar(x,par);
> if(prob > 0) logL += TMath::Log(prob) * weight[i];
> else logL += logEpsilon * weight[i];
> //////////////////////////////////////////////////////////
>
>
> Best Regards,
> Hajime
>
> > Hi Hajime,
> >
> > Thanks very much for your correction when using the weights.
> > iIt clearly improves the fit.
> >
> > It is now in CVS.
> >
> > Rene Brun
> >
> > Hajime Nanjyo wrote:
> > >
> > > Dear ROOTers,
> > >
> > > I found a bug in the function TreeUnbinnedFitLikelihood in TTreePlayer.cxx,
> > > which was related to a bug in TTree::UnbinnedFit.
> > >
> > > In the function TreeUnbinnedFitLikelihood, event weights were treated as
> > > follows.
> > > prob = fitfunc->EvalPar(x,par) * weight[i]/sum;
> > >
> > > I think it should be
> > > prob = TMath::Power(fitfunc->EvalPar(x,par),weight[i]/sum).
> > >
> > > So I propose to change the original code,
> > >
> > > prob = fitfunc->EvalPar(x,par) * weight[i]/sum;
> > > if(prob > 0) logL += TMath::Log(prob);
> > > else logL += logEpsilon;
> > >
> > > to
> > >
> > > prob = fitfunc->EvalPar(x,par);
> > > if(prob > 0) logL += TMath::Log(prob) * weight[i]/sum;
> > > else logL += logEpsilon * weight[i]/sum;
> > > .
> > >
> > > I show a small sample, "test.C" to indicate the problem below.
> > >
> > > It can be compiled and executed as follows.
> > > g++ -o test `root-config --cflags` `root-config --libs` -lMinuit test.C
> > > ./test
> > >
> > > In this code, 20000 events with weights are generated and fitted with two ways.
> > > 1. Fill a histogram with the weights and fit it.
> > > 2. TTree::UnbinnedFit with weights.
> > >
> > > The results should be same but not the same currently.
> > > I already checked that the fix of the TreeUnbinnedFitLikelihood
> > > remedy the problem.
> > >
> > > The ROOT under /afs/cern.ch/sw/root/v4.00.08/rh73_gcc32/root was used
> > > with gcc version 3.2 in Red Hat Linux release 7.3 (Valhalla)
> > >
> > > Could you treat the problem please.
> > >
> > > Best Regards,
> > > Hajime
> > >
> > > test.C
> > > //////////////////////////////////////////////////////////////
> > > #include <iostream>
> > > #include <cstdio>
> > > #include "TROOT.h"
> > > #include "TApplication.h"
> > > #include "TFile.h"
> > > #include "TTree.h"
> > > #include "TCanvas.h"
> > > #include "TF1.h"
> > > #include "TH1.h"
> > >
> > > int main(int argc,char** argv)
> > > {
> > > TApplication theApp("App", &argc, argv);
> > > gApplication->Init();
> > >
> > > //---------------------------------
> > > // 1.) generate n0 events following 1+x^2+8/3*a0*x with weight w0.
> > > // 2.) generate n1 events following 1+x^2+8/3*a1*x with weight w1.
> > > // 3.) fill (n0+n1) events in TTree tr0 in TFile test.root.
> > > // 4.) Fill values to the TH1D his with weight.
> > > // 5.) Fit the histogram with [0]*(1+x^2+8/3*[1]*x)
> > > // 6.) Performe UnbinnedFit with p.d.f of 1+x^2+8/3*[0]*x
> > > //---------------------------------
> > >
> > > int n0=10000;
> > > int n1=10000;
> > >
> > > double a0=0.3;
> > > double a1=0.6;
> > > double w0=0.4;
> > > double w1=0.8;
> > > //---------------------------------
> > >
> > >
> > > TFile tf("test.root","RECREATE");
> > > TF1 f0("f0","1+x*x+8./3.*[0]*x",-1,1);
> > > f0.SetParameter(0,a0);
> > > f0.SetNpx(1000);
> > > TF1 f1("f1","1+x*x+8./3.*[0]*x",-1,1);
> > > f1.SetParameter(0,a1);
> > > f1.SetNpx(1000);
> > >
> > > Double_t value, weight;
> > >
> > > TTree tr0("tr0","tr0");
> > > tr0.Branch("value",&value,"value/D");
> > > tr0.Branch("weight",&weight,"weight/D");
> > >
> > > for(int i=0;i<n0;i++) {
> > > value=f0.GetRandom();
> > > weight=w0;
> > > tr0.Fill();
> > > }
> > > for(int i=0;i<n1;i++) {
> > > value=f1.GetRandom();
> > > weight=w1;
> > > tr0.Fill();
> > > }
> > > tr0.Write();
> > > tf.Close();
> > >
> > > //---------------------------------
> > >
> > > TFile tfi("test.root");
> > > TTree* tr;
> > > tr = (TTree*)tfi.Get("tr0");
> > >
> > > TF1 fun("fun","[0]*(1+x*x+8./3.*[1]*x)",-1,1);
> > >
> > > TH1D his("his","his",40,-1,1);
> > > tr->Draw("value>>his","weight");
> > > his.Fit("fun","","",-1,1);
> > > Double_t A0= fun.GetParameter(1);
> > > //his.Draw("E");
> > > //gPad->Modified();
> > > //gPad->Update();
> > > //std::cout << "Hit Return" << std::endl;
> > > //std::getchar();
> > > //---------------------------------
> > >
> > > TF1 func("func","1+x*x+8./3.*[0]*x",-1,1);
> > > tr->UnbinnedFit("func","value","weight","VEM");
> > > Double_t A1= func.GetParameter(0);
> > >
> > > tfi.Close();
> > >
> > > std::cout << A0 << "\t" << A1 << std::endl;
> > > return 0;
> > > }
> > > //////////////////////////////////////////////////////////////
>
This archive was generated by hypermail 2b29 : Sun Jan 02 2005 - 05:50:09 MET