I have a (mostly) linear funtion to fit. Please help

demiller · August 27, 2005, 1:02am

Hello all.

I am fitting to a function of the following form:

f(x) = P[0] + P[1] * (x - P[3]) + P[2] * log(x - P[3])

or,

g(x) = P[0] + P[1] * (x) + P[2] * log(x)

h(x) = x - P[3]

So that,

f(x) = g(h(x))

So, you can see that it is linear in the parameters P[0,1 & 2],
but non-linear in parameter P[3].

I also see that we have a fast new TLinearFitter class.

What is the best way to nest the fitting routines for this
problem?

I imagine something like this:

I make a LinearFitter call for f(x) fitting P[3]
where the parameters for P[0-2] are LinearFitter
calculated for each trial value of P[3]. This way
the P[3] fitting call can treat the values P[0-2]
as constants and then simply do a 1-D fit for P[3].

This should give me a significant speed improvement.

So, in the new TFormula syntax,

g(x) = "1 ++ x ++ log(x) "

right? And a Fit call with this g(x) should be fast.

But,

f(x) = g(h(x)) = g(x-P[3]) is tricker.

I would like to use Root based fitting routines as much
as possible here.

Can anyone suggest a good way to do this?

Thanks.

-dave miller

anna · August 29, 2005, 10:57am

Hello, Dave,

I’m not sure I understand what you want to do. Is it: you have a good guess about the value of P[3], you want then to define g(x) as you wrote, make a linear fit and use the fitted parameters as initial values of P[0], P[1] and P[2] for Minuit fit? If so, you’ll have to define 2 functions, g(x) - the linear one 1++x++log(x), and a regular TF1 [0]+[1]*(x-[3])+[2]*log(x-[3]), in which you’d set the parameters according to the results of the linear fit. You can’t define a function as g(x-[3]).

Anna

demiller · October 20, 2005, 12:06am

Hello again.

What I need to do is a linear (in parameters) fit on C0-C2 where I control
the value of T0 externally (say to 2.5 for example).

The function is this: f(x)=C0+C1(x-T0)+C2(log(x-T0)).

So, I need to do a linear fit on
the function: f(x)=C0+C1(x-2.5)+C2(log(x-2.5)).

int main(void){

double T0=2.5; <- This is my externally controlled T0 value.

TVectorD *testVector=new TVectorD(1,2,3,4,5,“END”);

TH1D *testHisto=new TH1D(testVector);

TF1 outerTF1=new TF1(“outerTF1”,“1++(x-T0)++(log(x-T0)”,-2,20);

*testHisto->Fit(“outerTF1”);

}

But Root complains that T0 is not an arithemtic expression.

What am I doing wrong here?

brun · October 20, 2005, 5:06am

Replace T0 by [0]
and give the initial (or only value) of your parameter with
outerTF1->SetParameter(0,T0);
see examples in the fitting tutorials

Rene

anna · October 20, 2005, 10:48am

For the linear fit, you’ll have to put it as a constant, 1++(x-2.5)++… For now, the linear fitter can’t handle the case when there are other, nonlinear, parameters, even when they don’t have to be fit. I hope to add this functionality in the near future.

demiller · October 20, 2005, 1:46pm

Hello.

I can’t use a literal constant in the TFormula such as " 1++(x-2.5)++…"
because I will need to use yet to be determined values (not usually 2.5).

TFormula can accept c style functions.

Can I use a wrapper of some kind?

Such as,

namespace wrapper{
double internalVariable;
}
inline double variableWrapper(double*,double*){
reurn wrapper::internalVariable;
}

And then use the function,
TF1 *innerFit =
new TF1(“innerFit”,“1++(x-variableWrapper)++(log(x-variableWrapper))”,-2,20);

*testHisto->Fit(“innerFit”);

Or a TF1 wrapper around the c function?
Like this,
TF1 *varWrap = new TF1 (“varWrap”,variableWrapper,-2,20);
TF1 *innerFit =
new TF1 (“innerFit”,“1++(x-varWrap)++(log(x-varWrap))”,
-2,20)
*testHisto->Fit(“innerFit”);

Will anything like this work?

Thanks you.

demiller · October 20, 2005, 4:02pm

Here is another question.

Can I convert my double to a string. And then run the fit.

e.g.

double outerParameterDbl=123.45; // This is my externally set parameter.
stringstream outParStrstr;
outParStrstr << outerParameterDbl; // insert double into stream

string outParString;
outParStrstr >> outParString;
outParString = “(1++(x-” + outParString + “)++(log(x-” + outParString + “)))” ;

This string should now read: “(1++(x-123.45)++(log(x-123.45)))”

So,
TF1 *myTF1=new TF1(“myTF1”,outParString,-2,20);

And I call:
*myTH1D->Fit(“myTF1”);

Will this work?

anna · October 20, 2005, 4:29pm

Nice idea! Indeed, if you print a double into a string and then use this string as a title, it works.

demiller · November 4, 2005, 5:02pm

I am now having trouble invoking the linear fitter.

If I run the following code in cint I get the expected result:
$ root

```
                                    *
```
```
   W E L C O M E  to  R O O T       *
```
```
                                    *
```
Version 5.06/00 30 October 2005 *
```
                                    *
```
You are welcome to visit our Web site *

     [root.cern.ch](http://root.cern.ch)            *

```
                                    *
```

FreeType Engine v2.1.9 used to render TrueType fonts.
Compiled on 1 November 2005 for linux with thread support.

CINT/ROOT C/C++ Interpreter version 5.16.2, July 14 2005
Type ? for help. Commands must be C++ statements.
Enclose multiple statements between { }.

root [0] double dataArray[10]={1,2,3,4,5,6,7,8,9,10};

root [1] TH1D *h1LogGamma=new TH1D(“h1logGamma”,“Log Gamma Histogram;Bucket Number;ADC Log Value”,10,0.,11.);

root [2] h1LogGamma->FillN(10,dataArray,dataArray);

root [3] TF1 *innerGammaFit=new TF1(“innerGammaFit”,“1++x++log(x)”,0,20);

root [4] innerGammaFit->Print();

    innerGammaFit : 1++x++log(x) Ndim= 1, Npar= 3, Noper= 12

fExpr[0] = [0] action = 140 action param = 0
fExpr[1] = 1 action = 141 action param = 0
fExpr[2] = * action = 3 action param = 0
fExpr[3] = [1] action = 140 action param = 1
fExpr[4] = x action = 144 action param = 0
fExpr[5] = * action = 3 action param = 0
fExpr[6] = + action = 1 action param = 0
fExpr[7] = [2] action = 140 action param = 2
fExpr[8] = x action = 144 action param = 0
fExpr[9] = log action = 30 action param = 0
fExpr[10] = * action = 3 action param = 0
fExpr[11] = + action = 1 action param = 0
Optimized expresion
fExpr[0] = [0],1,* action = 148 action param = 0
fExpr[1] = [1],x,* action = 148 action param = 0
fExpr[2] = + action = 1 action param = 0
fExpr[3] = [2] action = 146 action param = 2
fExpr[4] = x,log action = 147 action param = 0
fExpr[5] = * action = 3 action param = 0
fExpr[6] = + action = 1 action param = 0
Par 0 p0 = 0
Par 1 p1 = 0
Par 2 p2 = 0

root [5] h1LogGamma->Fit(innerGammaFit);

Fitting results:
Parameters:
NO. VALUE ERROR
0 0.500000 1.034926
1 0.909091 0.526501
2 0.000000 1.475368
TCanvas::MakeDefCanvas: created default TCanvas with name c1

But when I compile the following code (testFit.cc):

#include
#include “TH1D.h”
#include “TF1.h”

int main(){

double dataArray[10]={1,2,3,4,5,6,7,8,9,10};

TH1D *h1LogGamma=new TH1D(“h1logGamma”,“Log Gamma Histogram;Bucket Number;ADC Log Value”,10,0.,11.);

h1LogGamma->FillN(10,dataArray,dataArray);

TF1 *innerGammaFit=new TF1(“innerGammaFit”,“1++x++log(x)”,0,20);

innerGammaFit->Print();

std::cout<<std::endl<<“Will this fit?”<<std::endl;
std::cin.get();

h1LogGamma->Fit(“innerGammaFit”);

std::cout<<std::endl<<“Did this fit?”<<std::endl;
std::cin.get();
}

With g++:
g++ -ggdb -o testFit root-config --cflags root-config --glibs testFit.cc

Where:
root-config --cflags = "-pthread -I/home/demiller/ROOT5.06.00/root/include"
root-config --glibs = “-L/home/demiller/ROOT5.06.00/root/lib -lCore -lCint -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lGui -pthread -lm -ldl -rdynamic”

I get a crash when I call the fit:

OUTPUT:
$ ./testFit
innerGammaFit : 1++x++log(x) Ndim= 1, Npar= 3, Noper= 12
fExpr[0] = [0] action = 140 action param = 0
fExpr[1] = 1 action = 141 action param = 0
fExpr[2] = * action = 3 action param = 0
fExpr[3] = [1] action = 140 action param = 1
fExpr[4] = x action = 144 action param = 0
fExpr[5] = * action = 3 action param = 0
fExpr[6] = + action = 1 action param = 0
fExpr[7] = [2] action = 140 action param = 2
fExpr[8] = x action = 144 action param = 0
fExpr[9] = log action = 30 action param = 0
fExpr[10] = * action = 3 action param = 0
fExpr[11] = + action = 1 action param = 0
Optimized expresion
fExpr[0] = [0],1,* action = 148 action param = 0
fExpr[1] = [1],x,* action = 148 action param = 0
fExpr[2] = + action = 1 action param = 0
fExpr[3] = [2] action = 146 action param = 2
fExpr[4] = x,log action = 147 action param = 0
fExpr[5] = * action = 3 action param = 0
fExpr[6] = + action = 1 action param = 0
Par 0 p0 = 0
Par 1 p1 = 0
Par 2 p2 = 0

Will this fit?

*** Break *** segmentation violation
Generating stack trace…
0x001ad761 in TH1::Fit(TF1*, char const*, char const*, double, double) + 0x509 from /home/demiller/ROOT5.06.00/root/lib/
libHist.so
0x001ad24f in TH1::Fit(char const*, char const*, char const*, double, double) + 0x1a7 from /home/demiller/ROOT5.06.00/root/
lib/libHist.so
0x08048dc2 in main + 0x1e2 from ./testFit
0x0427379d in __libc_start_main + 0xed from /lib/tls/libc.so.6
0x08048b55 in std::ios_base::Init::~Init() + 0x31 from ./testFit
Abort

Can anyone explain this to me?

Thank you,

-David Eric Miller

brun · November 4, 2005, 5:08pm

Could you send the shortest RUNNING script/file reproducing the problem?

Rene

demiller · November 4, 2005, 5:52pm

Here is code which shows the problem.

File: testFit.cc:

#include
#include “TH1D.h”
#include “TF1.h”

int main(){

double dataArray[10]={1,2,3,4,5,6,7,8,9,10};

TH1D *h1LogGamma=new TH1D(“h1logGamma”,“Log Gamma Histogram;Bucket Number;ADC Log Value”,10,0.,11.);

h1LogGamma->FillN(10,dataArray,dataArray);

TF1 *innerGammaFit=new TF1(“innerGammaFit”,“1++x++log(x)”,0,20);

innerGammaFit->Print();

std::cout<<std::endl<<“Will this fit?”<<std::endl;
std::cin.get();

h1LogGamma->Fit(“innerGammaFit”);

std::cout<<std::endl<<“Did this fit?”<<std::endl;
std::cin.get();
}

Compiled with:
g++ (GCC) 3.2.3 20030502 (Red Hat Linux 3.2.3-42)

g++ -ggdb -o testFit -pthread -I/home/demiller/ROOT5.06.00/root/include -L/home/demiller/ROOT5.06.00/root/lib -lCore -lCint -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lGui -pthread -lm -ldl -rdynamic testFit.cc

And the output from the executable is:

$ ./testFit
innerGammaFit : 1++x++log(x) Ndim= 1, Npar= 3, Noper= 12
fExpr[0] = [0] action = 140 action param = 0
fExpr[1] = 1 action = 141 action param = 0
fExpr[2] = * action = 3 action param = 0
fExpr[3] = [1] action = 140 action param = 1
fExpr[4] = x action = 144 action param = 0
fExpr[5] = * action = 3 action param = 0
fExpr[6] = + action = 1 action param = 0
fExpr[7] = [2] action = 140 action param = 2
fExpr[8] = x action = 144 action param = 0
fExpr[9] = log action = 30 action param = 0
fExpr[10] = * action = 3 action param = 0
fExpr[11] = + action = 1 action param = 0
Optimized expresion
fExpr[0] = [0],1,* action = 148 action param = 0
fExpr[1] = [1],x,* action = 148 action param = 0
fExpr[2] = + action = 1 action param = 0
fExpr[3] = [2] action = 146 action param = 2
fExpr[4] = x,log action = 147 action param = 0
fExpr[5] = * action = 3 action param = 0
fExpr[6] = + action = 1 action param = 0
Par 0 p0 = 0
Par 1 p1 = 0
Par 2 p2 = 0

Will this fit?

*** Break *** segmentation violation
Generating stack trace…
0x001ad761 in TH1::Fit(TF1*, char const*, char const*, double, double) + 0x509 from /home/demiller/ROOT5.06.00/root/lib/
libHist.so
0x001ad24f in TH1::Fit(char const*, char const*, char const*, double, double) + 0x1a7 from /home/demiller/ROOT5.06.00/root/
lib/libHist.so
0x08048dc2 in main + 0x1e2 from ./testFit
0x0427379d in __libc_start_main + 0xed from /lib/tls/libc.so.6
0x08048b55 in std::ios_base::Init::~Init() + 0x31 from ./testFit
Abort

Thank you.

-David Eric Miller

demiller · November 4, 2005, 5:58pm

Attached is the executeable from the code in my previous post.

demiller · November 4, 2005, 6:00pm

It looks like the attachemnt failed.

Hopefully you can compile this test code.

brun · November 4, 2005, 6:32pm

You did not attach the file(s).
Types must be x.C, x.tar.gz

Rene

demiller · November 4, 2005, 6:33pm

Here is some output from gdb:

Program received signal SIGSEGV, Segmentation fault.
0x0045fb4d in TClass::New () from /home/demiller/ROOT5.06.00/root/lib/libCore.so

(gdb) where

#0 0x0045fb4d in TClass::New () from /home/demiller/ROOT5.06.00/root/lib/libCore.so
#1 0x00a29761 in TH1::Fit () from /home/demiller/ROOT5.06.00/root/lib/libHist.so
#2 0x08048dc0 in main () at testFit.cc:20

-David Miller

demiller · November 4, 2005, 6:38pm

Here is the code.
testFit.gz (10 KB)
testFit.cc.gz (321 Bytes)

brun · November 4, 2005, 7:51pm

I do not see a problem with your trivial test. It produces the output in the attachement.
In your program, I suggest to replace
int main() {
by
int testFit() {

and from a ROOT interactive session, do
root > .x testFit.C

where testFit.C is the modified test program.

Rene

demiller · November 4, 2005, 9:36pm

Hello Rene,

An interactive ROOT session has always executed the macro correctly.

It is the COMPILED code which fails to run properly.

Did you execute the program I sent you?

It seg faults inside the Fit call.

Can you COMPILE my test code and then run it without
the segmentation fault?

Again, here is the gdb output from the COMPILED program:

Program received signal SIGSEGV, Segmentation fault.
0x0045fb4d in TClass::New () from /home/demiller/ROOT5.06.00/root/lib/libCore.so

(gdb) where

#0 0x0045fb4d in TClass::New () from /home/demiller/ROOT5.06.00/root/lib/libCore.so
#1 0x00a29761 in TH1::Fit () from /home/demiller/ROOT5.06.00/root/lib/libHist.so
#2 0x08048dc0 in main () at testFit.cc:20

Again,
Thank you verry much.

-David Miller

demiller · November 4, 2005, 10:03pm

Here is some more output from my test program from the line
h1LogGamma->Fit(innerGammaFit);

*** Break *** segmentation violation
Generating stack trace…
0x00501140 in TH1::Fit(TF1*, char const*, char const*, double, double) + 0x36e from /usr/local/root/pro/lib/libHist.so
0x08048dc4 in main + 0x1e0 from ./runTest
0x0743e78a in __libc_start_main + 0xda from /lib/tls/libc.so.6
0x08048b55 in std::ios_base::Init::~Init in-charge + 0x31 from ./runTest
Abort (core dumped)

And attached is the core file.

-David Miller
core.3928.gz (602 KB)

brun · November 4, 2005, 10:14pm

When you run with your own main program, you must run with a TApplication object instantiated. ie:

[code]#include
#include “TApplication.h”
#include “TH1D.h”
#include “TF1.h”

int main(int argc, char **argv)
{
TApplication theApp(“App”, &argc, argv);

double dataArray[10]={1,2,3,4,5,6,7,8,9,10};

TH1D *h1LogGamma=new TH1D(“h1logGamma”,“Log Gamma Histogram;Bucket Number;ADC Log Value”,10,0.,11.);

h1LogGamma->FillN(10,dataArray,dataArray);

TF1 *innerGammaFit=new TF1(“innerGammaFit”,“1++x++log(x)”,0,20);

innerGammaFit->Print();

std::cout<<std::endl<<“Will this fit?”<<std::endl;
std::cin.get();

h1LogGamma->Fit(innerGammaFit);

std::cout<<std::endl<<“Did this fit?”<<std::endl;
std::cin.get();
return 0;
}[/code]

Rene