Go to the documentation of this file.
1 /**********************************************************************
2  * File: quadlsq.cpp (Formerly qlsq.c)
3  * Description: Code for least squares approximation of quadratics.
4  * Author: Ray Smith
5  * Created: Wed Oct 6 15:14:23 BST 1993
6  *
7  * (C) Copyright 1993, Hewlett-Packard Ltd.
9  ** you may not use this file except in compliance with the License.
10  ** You may obtain a copy of the License at
12  ** Unless required by applicable law or agreed to in writing, software
14  ** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15  ** See the License for the specific language governing permissions and
16  ** limitations under the License.
17  *
18  **********************************************************************/
19
20 #include <stdio.h>
21 #include <math.h>
23 #include "tprintf.h"
24
25 // Minimum variance in least squares before backing off to a lower degree.
26 const double kMinVariance = 1.0 / 1024;
27
28 /**********************************************************************
29  * QLSQ::clear
30  *
31  * Function to initialize a QLSQ.
32  **********************************************************************/
33
34 void QLSQ::clear() { // initialize
35  a = 0.0;
36  b = 0.0;
37  c = 0.0;
38  n = 0; // No elements.
39  sigx = 0.0; // Zero accumulators.
40  sigy = 0.0;
41  sigxx = 0.0;
42  sigxy = 0.0;
43  sigyy = 0.0;
44  sigxxx = 0.0;
45  sigxxy = 0.0;
46  sigxxxx = 0.0;
47 }
48
49
50 /**********************************************************************
52  *
53  * Add an element to the accumulator.
54  **********************************************************************/
55
56 void QLSQ::add(double x, double y) {
57  n++; // Count elements.
58  sigx += x; // Update accumulators.
59  sigy += y;
60  sigxx += x * x;
61  sigxy += x * y;
62  sigyy += y * y;
63  sigxxx += static_cast<long double>(x) * x * x;
64  sigxxy += static_cast<long double>(x) * x * y;
65  sigxxxx += static_cast<long double>(x) * x * x * x;
66 }
67
68
69 /**********************************************************************
70  * QLSQ::remove
71  *
72  * Delete an element from the accumulator.
73  **********************************************************************/
74
75 void QLSQ::remove(double x, double y) {
76  if (n <= 0) {
77  tprintf("Can't remove an element from an empty QLSQ accumulator!\n");
78  return;
79  }
80  n--; // Count elements.
81  sigx -= x; // Update accumulators.
82  sigy -= y;
83  sigxx -= x * x;
84  sigxy -= x * y;
85  sigyy -= y * y;
86  sigxxx -= static_cast<long double>(x) * x * x;
87  sigxxy -= static_cast<long double>(x) * x * y;
88  sigxxxx -= static_cast<long double>(x) * x * x * x;
89 }
90
91
92 /**********************************************************************
93  * QLSQ::fit
94  *
95  * Fit the given degree of polynomial and store the result.
96  * This creates a quadratic of the form axx + bx + c, but limited to
97  * the given degree.
98  **********************************************************************/
99
100 void QLSQ::fit(int degree) {
101  long double x_variance = static_cast<long double>(sigxx) * n -
102  static_cast<long double>(sigx) * sigx;
103
104  // Note: for computational efficiency, we do not normalize the variance,
105  // covariance and cube variance here as they are in the same order in both
106  // nominators and denominators. However, we need be careful in value range
107  // check.
108  if (x_variance < kMinVariance * n * n || degree < 1 || n < 2) {
109  // We cannot calculate b reliably so forget a and b, and just work on c.
110  a = b = 0.0;
111  if (n >= 1 && degree >= 0) {
112  c = sigy / n;
113  } else {
114  c = 0.0;
115  }
116  return;
117  }
118  long double top96 = 0.0; // Accurate top.
119  long double bottom96 = 0.0; // Accurate bottom.
120  long double cubevar = sigxxx * n - static_cast<long double>(sigxx) * sigx;
121  long double covariance = static_cast<long double>(sigxy) * n -
122  static_cast<long double>(sigx) * sigy;
123
124  if (n >= 4 && degree >= 2) {
125  top96 = cubevar * covariance;
126  top96 += x_variance * (static_cast<long double>(sigxx) * sigy - sigxxy * n);
127
128  bottom96 = cubevar * cubevar;
129  bottom96 -= x_variance *
130  (sigxxxx * n - static_cast<long double>(sigxx) * sigxx);
131  }
132  if (bottom96 >= kMinVariance * n * n * n * n) {
133  // Denominators looking good
134  a = top96 / bottom96;
135  top96 = covariance - cubevar * a;
136  b = top96 / x_variance;
137  } else {
138  // Forget a, and concentrate on b.
139  a = 0.0;
140  b = covariance / x_variance;
141  }
142  c = (sigy - a * sigxx - b * sigx) / n;
143 }
void remove(double x, double y)