4 template <
class T>
T sqr(
T t) {
return t*
t;}
22 if (hasErrors)
return theResult;
23 if (hasValues && !doErrors)
return theResult;
25 double F0, F1, F2, F3, F4, F0y, F1y, F2y;
26 F0 = F1 = F2 = F3 = F4 = F0y = F1y = F2y = 0.;
28 typedef vector<Point>::const_iterator
IT;
37 F1 += w*x; F1y += w*x*y;
38 pow = x*x; F2 += w*
pow; F2y += w*pow*y;
39 pow *= x; F3 += w*
pow;
40 pow *= x; F4 += w*
pow;
43 Column cA = { F0, F1, F2 };
44 Column cB = { F1, F2, F3 };
45 Column cC = { F2, F3, F4 };
46 Column cY = { F0y, F1y, F2y };
48 double det0 = det(cA, cB, cC);
51 theResult.parA = det(cY, cB, cC) / det0;
52 theResult.parB = det(cA, cY, cC) / det0;
53 theResult.parC = det(cA, cB, cY) / det0;
55 Column cCY = {F0y-theResult.parC*F2, F1y-theResult.parC*F3, F2y-theResult.parC*F4};
56 double det0C = det(cA, cB);
57 theResult.parA = det(cCY, cB) / det0C;
58 theResult.parB = det(cA, cCY) / det0C;
62 double vAA, vBB, vCC, vAB, vAC, vBC;
63 vAA = vBB = vCC = vAB = vAC = vBC = 0.;
66 if (!doErrors)
return theResult;
68 if( !hasWeights &&
dof() > 0) {
70 double scale_w = 1./
chi2()/
dof();
71 for (IT ip =
points.begin(); ip !=
points.end(); ip++) ip->w *= scale_w;
80 double dXBC = det(cX, cB, cC);
81 double dAXC = det(cA, cX, cC);
82 double dABX = det(cA, cB, cX);
87 vAB += w * dXBC * dAXC;
88 vAC += w * dXBC * dABX;
89 vBC += w * dAXC * dABX;
92 theResult.varAA = vAA/
sqr(det0);
93 theResult.varBB = vBB/
sqr(det0);
94 theResult.varCC = vCC/
sqr(det0);
95 theResult.varAB = vAB/
sqr(det0);
96 theResult.varAC = vAC/
sqr(det0);
97 theResult.varBC = vBC/
sqr(det0);
107 for ( vector<Point>::const_iterator
109 mychi2 += ip->w *
sqr(ip->y - fun(ip->x));
116 return theResult.parA + theResult.parB*x + theResult.parC*x*x;
121 int nPar = 3;
if (hasFixedParC) nPar--;
123 return (nDof > 0) ? nDof : 0;
129 return c1.
r1 * c2.
r2 * c3.
r3
const Result & result(bool doErrors) const
double det(const Column &c1, const Column &c2, const Column &c3) const
void addPoint(double x, double y)
std::vector< LinkConnSpec >::const_iterator IT
double fun(double x) const
Power< A, B >::type pow(const A &a, const B &b)