Linear Fitting @ >3rd order

Hi,

With the macro below, I find that having the line with 4th order terms in my TF2 causes it to malfunction during a fit. While my input paramters are all in the range -5 => 5, the fitter converges outside this range with a very high chi-square.

If I comment out the line with 4th order terms, everything looks ok. Is this a known limitation with the linear fit algorithm at high order (i.e singular matrices)?

In addition, when I add in the 5th order line, I get the following error:

Error in TFormula::Compile: ‘)’ is expected
Error in TF2::TF2: function: func/1 ++ x ++ y ++ x^2 ++ x y ++ y^2 ++ x^3 ++ x^2y ++ x y^2 ++ y^3 ++ x^4 ++ x^3y ++ x^2y^2 ++ x y^3 ++ y^4++ x^5 ++ x^4y ++ x^3y^2 ++ x^2*y^3 ++ x *y^4 ++ y^5 has 0 parameters instead of 2

Thanks for any help.

Cheers,
Chris[code]void testFit() {

const Int_t n(100);

Double_t min(0),max(5000);
Double_t x[n],y[n],z[n];

TF2 func(“func”,
“1 “
”++ x ++ y “
”++ x^2 ++ x y ++ y^2 “
”++ x^3 ++ x^2
y ++ x y^2 ++ y^3 “
”++ x^4 ++ x^3
y ++ x^2y^2 ++ x y^3 ++ y^4"
//"++ x^5 ++ x^4
y ++ x^3
y^2 ++ x^2*y^3 ++ x *y^4 ++ y^5”
, min, max, min, max);

TRandom rand;

cout << "Input params " << endl;
for (Int_t i=0; i<func.GetNpar(); ++i ) {
Double_t value=rand.Uniform(-5,5);
cout << i << " " << value << endl;
func.SetParameter( i, value );
}

for (Int_t i=0; i<n; ++i ) {
x[i]= rand.Uniform(min,max);
y[i]= rand.Uniform(min,max);
z[i]= func.Eval(x[i],y[i]) + rand.Gaus(0,1);
}

TGraph2D graph(n,x,y,z);

graph.Fit(&func, “”);

// print the output – errors are zero w/ linear fitter :-/
for ( Int_t i=0; i<func.GetNpar(); ++i ) {
printf(“par[%2i] %9.4g %9.4g\n”, i,
func.GetParameter(i),
func.GetParError(i));
}
// dump the chisquare, too
printf(“chisquare= %9.4f\n”, func.GetChisquare() );

}[/code]

Hi Chris,

If you try, for example, to center your data and move its range from (0, 5000) to (-2500, 2500), the results will be much better. But in general, of course, singular design matrices are a limitation of the linear fitter algorithm, as they are for most fitting algorithms. There are ways to go around it, but it depends on the data and on the fitting function.
As for your 5th order function, it doesn’t compile, because after some linear function processing its length exceeds the internal TFormula limit of 256 symbols. The only thing I can recommend you in this case is to use the linear fitter directly, where there is no such limitation.

Hi Anna,

Thanks for your response. [By the way, I’m using ROOT 5.08/00 on Debian linux).

I’m glad you could verify it was singular matrices that were causing the problem. Do you think it would be useful if there was a warning message printed by the fitter in these cases, to warn the user?

I switched to the linear fitter now, so I generate data with a TF2 and fit it with a TLinearFitter - a little awkward, but it seems to work. The problem is that I get a crash with EvalRobust, as I reported in another thread. See the traceback & macro below.

I was also curious what if there is a way to get parameter errors and chisquare when using the linear fitter. Using the macro below with “Eval” instead of “EvalRobust” works great, but the algorithm doesn’t compute errors directly.

Cheers,
Chris

root [6] .x linearFit.C
Input params
0 -3.92085
1 0.874751
2 -1.81027
3 -3.53562
4 -2.00777
5 -4.59299

*** Break *** segmentation violation
Generating stack trace…
0xffffe420 in
0x40b5918a in TFormula::EvalPar(double const*, double const*) + 0x4e from /usr/lib/root/libHist.so.5.08
0x41fa7413 in TLinearFitter::CStep(int, int, double*, int*, int*, int, int) + 0xe2d from /usr/lib/root/5.08/libMinuit.so
0x41fa49d5 in TLinearFitter::EvalRobust(double) + 0x281 from /usr/lib/root/5.08/libMinuit.so
0x41fce154 in from /usr/lib/root/5.08/libMinuit.so
0x40880571 in G__ExceptionWrapper + 0x5d from /usr/lib/root/libCint.so.5.08
0x409339a4 in G__call_cppfunc + 0x30e from /usr/lib/root/libCint.so.5.08
0x409207eb in G__interpret_func + 0x8ab from /usr/lib/root/libCint.so.5.08
0x40906989 in G__getfunction + 0x1c3d from /usr/lib/root/libCint.so.5.08
0x409a5e74 in G__getstructmem + 0xb9e from /usr/lib/root/libCint.so.5.08
0x4099d0fc in G__getvariable + 0x7b8 from /usr/lib/root/libCint.so.5.08
0x408fbd0c in G__getitem + 0x6dc from /usr/lib/root/libCint.so.5.08
0x408fa22a in G__getexpr + 0xc27e from /usr/lib/root/libCint.so.5.08
0x4094e9e2 in G__exec_function + 0x202 from /usr/lib/root/libCint.so.5.08
0x40956aef in G__exec_statement + 0x2a47 from /usr/lib/root/libCint.so.5.08
0x409223bc in G__interpret_func + 0x247c from /usr/lib/root/libCint.so.5.08
0x409070aa in G__getfunction + 0x235e from /usr/lib/root/libCint.so.5.08
0x408fbd97 in G__getitem + 0x767 from /usr/lib/root/libCint.so.5.08
0x408fa22a in G__getexpr + 0xc27e from /usr/lib/root/libCint.so.5.08
0x408edad2 in G__calc_internal + 0x312 from /usr/lib/root/libCint.so.5.08
0x4095d19b in G__process_cmd + 0x2359 from /usr/lib/root/libCint.so.5.08
0x4021b9cc in TCint::ProcessLine(char const*, TInterpreter::EErrorCode*) + 0x136 from /usr/lib/root/libCore.so.5.08
0x4021bb5d in TCint::ProcessLineSynch(char const*, TInterpreter::EErrorCode*) + 0x4f from /usr/lib/root/libCore.so.5.08
0x40150319 in TApplication::ProcessFile(char const*, int*) + 0x9d5 from /usr/lib/root/libCore.so.5.08
0x4014f8c1 in TApplication::ProcessLine(char const*, bool, int*) + 0x5f9 from /usr/lib/root/libCore.so.5.08
0x414ae978 in TRint::HandleTermInput() + 0x24e from /usr/lib/root/libRint.so.5.08
0x414ad205 in TTermInputHandler::Notify() + 0x27 from /usr/lib/root/libRint.so.5.08
0x414af5aa in TTermInputHandler::ReadNotify() + 0x14 from /usr/lib/root/libRint.so.5.08
0x402bb852 in TUnixSystem::CheckDescriptors() + 0x166 from /usr/lib/root/libCore.so.5.08
0x402ba697 in TUnixSystem::DispatchOneEvent(bool) + 0x185 from /usr/lib/root/libCore.so.5.08
0x401d0f44 in TSystem::InnerLoop() + 0x22 from /usr/lib/root/libCore.so.5.08
0x401d0eda in TSystem::Run() + 0x80 from /usr/lib/root/libCore.so.5.08
0x4015061b in TApplication::Run(bool) + 0x37 from /usr/lib/root/libCore.so.5.08 0x414ae2de in TRint::Run(bool) + 0x44c from /usr/lib/root/libRint.so.5.08
0x08048e64 in main + 0x90 from /usr/bin/root.exe
0x4161de36 in __libc_start_main + 0xc6 from /lib/libc.so.6
0x08048d31 in TApplicationImp::ShowMembers(TMemberInspector&, char*) + 0x3d from /usr/bin/root.exe
Root > Function linearFit() busy flag cleared

void linearFit() {

const Char_t *funcString=
"1 “
”++ x ++ y “
”++ x^2 ++ x *y ++ y^2 ";

const Int_t n(100);

Double_t min(0),max(5000);

TF2 func(“func”,funcString,min,max,min,max);

TRandom rand;

cout << "Input params " << endl;
for (Int_t i=0; i<func.GetNpar(); ++i ) {
Double_t value=rand.Uniform(-5,5);
cout << i << " " << value << endl;
func.SetParameter( i, value );
}

TLinearFitter fitter(2,funcString);
Double_t x[2],y;
for (Int_t i=0; i<n; ++i ) {
x[0]= rand.Uniform(min,max);
x[1]= rand.Uniform(min,max);
y= func.Eval(x[0],x[1]) + rand.Gaus(0,10);
fitter.AddPoint(x,y);
}

fitter.EvalRobust();

// print the output – errors are zero w/ linear fitter :-/
for ( Int_t i=0; i<func.GetNpar(); ++i ) {
printf(“par[%2i] %9.4g %9.4g\n”, i,
func.GetParameter(i),
func.GetParError(i));
}
// dump the chisquare, too
printf(“chisquare= %9.4f\n”, func.GetChisquare() );

}

Hi,

I got 5th order generation/fitting working with the example below, using an outside function for the TF2. However, I get a segmentation violation when the script ends. The output and macro are below. I’m using ROOT 5.08/00 on Debian linux.

Chris

root [0] .x linearFit.C
Input params
0 -3.92085
1 0.874751
2 -1.81027
3 -3.53562
4 -2.00777
5 -4.59299
6 -3.47878
7 4.56471
8 -0.201092
9 0.815624
10 4.4044
11 -2.58997
12 3.39339
13 -1.91606
14 -0.293188
15 -0.141964
16 4.75405
17 -2.441
18 2.27815
19 -0.415282
20 -3.0927
par[ 0] -3.921 0
par[ 1] 0.8748 0
par[ 2] -1.81 0
par[ 3] -3.536 0
par[ 4] -2.008 0
par[ 5] -4.593 0
par[ 6] -3.479 0
par[ 7] 4.565 0
par[ 8] -0.2011 0
par[ 9] 0.8156 0
par[10] 4.404 0
par[11] -2.59 0
par[12] 3.393 0
par[13] -1.916 0
par[14] -0.2932 0
par[15] -0.142 0
par[16] 4.754 0
par[17] -2.441 0
par[18] 2.278 0
par[19] -0.4153 0
par[20] -3.093 0
chisquare= 0.0000

*** Break *** segmentation violation
Killed

Double_t funcIt(Double_t *dep, Double_t *par ) {
Double_t x=dep[0];
Double_t y=dep[1];

Double_t x2=xx;
Double_t x3=x2
x;
Double_t x4=x3x;
Double_t x5=x4
x;

Double_t y2=yy;
Double_t y3=y2
y;
Double_t y4=y3y;
Double_t y5=y4
y;

return
par[0]+
+ par[1]* x + par[2]* y
+ par[3]* x2 + par[4]* x y + par[5] y2
+ par[6]* x3 + par[7]* x2y + par[8] x y2 + par[9] y3
+ par[10]* x4 + par[11]* x3y + par[12] x2y2 + par[13] x y3 + par[14] y4
+ par[15]* x5 + par[16]* x4y + par[17] x3y2 + par[18] x2y3 + par[19] x y4 + par[20] y5;
}

void linearFit() {

const Char_t funcString=
"1 “
”++ x ++ y “
”++ x^2 ++ x y ++ y^2 “
”++ x^3 ++ x^2
y ++ x y^2 ++ y^3 “
”++ x^4 ++ x^3
y ++ x^2
y^2 ++ x y^3 ++ y^4"
"++ x^5 ++ x^4
y ++ x^3y^2 ++ x^2y^3 ++ x *y^4 ++ y^5";

const Int_t n(200);

Double_t min(-2500),max(2500);

TF2 func(“func”,&funcIt,min,max,min,max,21);

TRandom rand;

cout << "Input params " << endl;
for (Int_t i=0; i<func.GetNpar(); ++i ) {
Double_t value=rand.Uniform(-5,5);
cout << i << " " << value << endl;
func.SetParameter( i, value );
}

TLinearFitter fitter(2,funcString);
Double_t x[2],y;
for (Int_t i=0; i<n; ++i ) {
x[0]= rand.Uniform(min,max);
x[1]= rand.Uniform(min,max);
y= func.Eval(x[0],x[1]) + rand.Gaus(0,10);
fitter.AddPoint(x,y);
}

fitter.Eval();

// print the output – errors are zero w/ linear fitter :-/
for ( Int_t i=0; i<func.GetNpar(); ++i ) {
printf(“par[%2i] %9.4g %9.4g\n”, i,
func.GetParameter(i),
func.GetParError(i));
}
// dump the chisquare, too
printf(“chisquare= %9.4f\n”, func.GetChisquare() );

}

Chris,

The EvalRobust crash was a bug that I fixed, so if you take the cvs head version it shouldn’t be a problem any more. As stated somewhere in the documentation, the least trimmed squares robust fitting algorithm doesn’t compute the errors on parameters, so I can’t help you here.
I’ll look into your last crash tomorrow

Cheers,
Anna

Hi Chris,

The cause of your second crash is fixed in the cvs head version, so everything should run smoothely now.
Please tell me, if you find something else not working :wink:

Anna