1 // leastsqs.c -- Implements a simple linear least squares best fit routine
3 // Written by Curtis Olson, started September 1997.
5 // Copyright (C) 1997 Curtis L. Olson - http://www.flightgear.org/~curt
7 // This library is free software; you can redistribute it and/or
8 // modify it under the terms of the GNU Library General Public
9 // License as published by the Free Software Foundation; either
10 // version 2 of the License, or (at your option) any later version.
12 // This library is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 // Library General Public License for more details.
17 // You should have received a copy of the GNU General Public License
18 // along with this program; if not, write to the Free Software
19 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
27 #include "leastsqs.hxx"
35 n*sum(xi*yi) - (sum(xi)*sum(yi))
36 b1 = --------------------------------
37 n*sum(xi^2) - (sum(xi))^2
40 b0 = sum(yi)/n - b1*(sum(xi)/n)
43 double sum_xi, sum_yi, sum_xi_2, sum_xi_yi;
46 void least_squares(double *x, double *y, int n, double *m, double *b) {
49 sum_xi = sum_yi = sum_xi_2 = sum_xi_yi = 0.0;
52 for ( i = 0; i < n; i++ ) {
55 sum_xi_2 += x[i] * x[i];
56 sum_xi_yi += x[i] * y[i];
59 /* printf("sum(xi)=%.2f sum(yi)=%.2f sum(xi^2)=%.2f sum(xi*yi)=%.2f\n",
60 sum_xi, sum_yi, sum_xi_2, sum_xi_yi); */
62 *m = ( (double)sum_n * sum_xi_yi - sum_xi * sum_yi ) /
63 ( (double)sum_n * sum_xi_2 - sum_xi * sum_xi );
64 *b = (sum_yi / (double)sum_n) - (*m) * (sum_xi / (double)sum_n);
66 /* printf("slope = %.2f intercept = %.2f\n", *m, *b); */
70 /* incrimentally update existing values with a new data point */
71 void least_squares_update(double x, double y, double *m, double *b) {
79 /* printf("sum(xi)=%.2f sum(yi)=%.2f sum(xi^2)=%.2f sum(xi*yi)=%.2f\n",
80 sum_xi, sum_yi, sum_xi_2, sum_xi_yi); */
82 *m = ( (double)sum_n * sum_xi_yi - sum_xi * sum_yi ) /
83 ( (double)sum_n * sum_xi_2 - sum_xi * sum_xi );
84 *b = (sum_yi / (double)sum_n) - (*m) * (sum_xi / (double)sum_n);
86 /* printf("slope = %.2f intercept = %.2f\n", *m, *b); */
91 return the least squares error:
97 double least_squares_error(double *x, double *y, int n, double m, double b) {
103 for ( i = 0; i < n; i++ ) {
104 error = y[i] - (m * x[i] + b);
105 sum += error * error;
106 // printf("%.2f %.2f\n", error, sum);
109 return ( sum / (double)n );
114 return the maximum least squares error:
118 double least_squares_max_error(double *x, double *y, int n, double m, double b){
120 double error, max_error;
124 for ( i = 0; i < n; i++ ) {
125 error = y[i] - (m * x[i] + b);
126 error = error * error;
127 if ( error > max_error ) {
132 return ( max_error );