Evaluate system of equations for the Moments Method

[marcogobbo]

Hello, I’m new to RooFit, trying to do various parameter estimation exercises involving minimization via LS (Binned), ML, EML and MM. While for the first two I used the additional parameter L (for likelihood) in the Fit function of TH1D, for the other two I used RooFit. For the EML, I wrote this code

    x_roo = RooRealVar("x_roo", "random variable", -1, 1)
    alpha_roo = RooRealVar("alpha_roo", "alpha", 0.40, 0.10, 0.70)
    beta_roo = RooRealVar("beta_roo", "beta", 0.60, 0.10, 0.70)
    scatter_pdf = RooGenericPdf("scatter_pdf", "scatter_pdf", "(1+alpha_roo*x_roo+beta_roo*(x_roo**2))/(2+2*beta_roo/3)", RooArgList(x_roo, alpha_roo, beta_roo))
    N_roo = RooRealVar("N_roo", "N_roo", 40000, 60000)
    ext_scatter_pdf = RooExtendPdf("ext_scatter_pdf", "extended scatter pdf", scatter_pdf, N_roo)

    hist = RooDataHist("data", "data histogram", RooArgList(x_roo), h3)
    r = ext_scatter_pdf.fitTo(hist, RooFit.Range("r1"))

    xframe = x_roo.frame(RooFit.Title("   "), RooFit.Range(-1.5, 1.5))
    hist.plotOn(xframe, RooFit.FillColor(kAzure + 7), RooFit.Binning(50))
    ext_scatter_pdf.plotOn(xframe, RooFit.LineColor(kRed + 1), RooFit.LineWidth(5))

    c3.cd()
    c3.SetTicks()
    c3.SetGrid()
    xframe.Draw()
    c3.Modified()
    c3.Update()
    c3.SaveAs("test_dist_3.pdf")

For the MM I don’t know how to proceed, or rather, I initialized the various elements, but now I should solve a system of equations of the following type (see image) from which to extract the alpha and beta values.

    x_roo = RooRealVar("x_roo", "random variable", -1, 1)
    alpha_roo = RooRealVar("alpha_roo", "alpha", 0.40, 0.10, 0.70)
    beta_roo = RooRealVar("beta_roo", "beta", 0.60, 0.10, 0.70)
    scatter_pdf = RooGenericPdf("scatter_pdf", "scatter_pdf", "(1+alpha_roo*x_roo+beta_roo*(x_roo**2))/(2+2*beta_roo/3)", RooArgList(x_roo, alpha_roo, beta_roo))

    integral_1 = RooRealIntegral("integral_1", "integral_1", "x_roo*scatter_pdf")
    integral_2 = RooRealIntegral("integral_2", "integral_2", "x_roo*x_roo*scatter_pdf")

How should I move? Thanks for your attention!

[eguiraud]

Hi @marcogobbo ,

we need @jonas 's help, let’s ping him.

Cheers,
Enrico

[jonas]

Hi @marcogobbo,

good luck with that! The method of moments is not used much in High Energy Physics, so it is not implemented in RooFit like maximum likelihood is.

However, at least you can get some objects to represent the moments relatively easily, no need to create the integrals yourself. To get the moments from the RooDataHist, you can use RooAbsData::moment(), which returns you double with the central moment you want to calculate.

For the pdf, you have RooAbsReal::moment(), which returns you a new RooFit object that represents a moment of given order.

Given these moments, it’s not clear to me how you want to proceed? Since you don’t have the moments analytically, how do you want to solve the system of equations? Again with least squares? You can only use RooFit for parameter estimate if you can write some test statistic that is then minimized. If you tell me what exactly you want to minimize then, I can maybe help you implement you this test statistic in RooFit.

Cheers,
Jonas

[marcogobbo]

The idea is actually quite simple from a mathematical point of view. My problem is understanding if I can solve this system of equations in a computational way, using RooFit functions, or if I am forced to evaluate the (simple) integral, add the random variables up to the number of events, and proceed through an error propagation.

double P = 0;
double Q = 0;

double NEV = 50000.;

for (int i=0; i<NEV; ++i){
P = P + distr_smearedV [i];
Q = Q + pow( distr_smearedV [i] ,2);
}

double errP = (1/ NEV )* sqrt (NEV *pow (0.05 ,2) );
double errQ = (1/ NEV ) *2*0.05* sqrt (Q);
P = P/NEV;
Q = Q/NEV;

double beta_hat = (2*Q -(2./3.) ) *(15/(6 -10* Q));
double beta_err = (20./ pow (3 -5*Q ,2) )* errQ ;
double alpha_hat = P*( beta_hat +3) ;
double alpha_err = sqrt (( pow ( beta_hat +3 ,2)*pow (errP ,2) )+( pow (P ,2) )*pow (
beta_err ,2) );

cout << alpha_hat << " +- " << alpha_err << endl ;
cout << beta_hat << " +- " << beta_err << endl ;

A friend gave me her analytical method (code above) if we want to call it that. At the same time, I was thinking if there was a way to implement this procedure with RooFit. If there is no immediate solution, I will also proceed in this way.

Thanks for the previous answer!