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 Library General Public
18 // License along with this library; if not, write to the
19 // Free Software Foundation, Inc., 59 Temple Place - Suite 330,
20 // Boston, MA 02111-1307, USA.
28 #include "leastsqs.hxx"
36 n*sum(xi*yi) - (sum(xi)*sum(yi))
37 b1 = --------------------------------
38 n*sum(xi^2) - (sum(xi))^2
41 b0 = sum(yi)/n - b1*(sum(xi)/n)
44 double sum_xi, sum_yi, sum_xi_2, sum_xi_yi;
47 void least_squares(double *x, double *y, int n, double *m, double *b) {
50 sum_xi = sum_yi = sum_xi_2 = sum_xi_yi = 0.0;
53 for ( i = 0; i < n; i++ ) {
56 sum_xi_2 += x[i] * x[i];
57 sum_xi_yi += x[i] * y[i];
60 /* printf("sum(xi)=%.2f sum(yi)=%.2f sum(xi^2)=%.2f sum(xi*yi)=%.2f\n",
61 sum_xi, sum_yi, sum_xi_2, sum_xi_yi); */
63 *m = ( (double)sum_n * sum_xi_yi - sum_xi * sum_yi ) /
64 ( (double)sum_n * sum_xi_2 - sum_xi * sum_xi );
65 *b = (sum_yi / (double)sum_n) - (*m) * (sum_xi / (double)sum_n);
67 /* printf("slope = %.2f intercept = %.2f\n", *m, *b); */
71 /* incrimentally update existing values with a new data point */
72 void least_squares_update(double x, double y, double *m, double *b) {
80 /* printf("sum(xi)=%.2f sum(yi)=%.2f sum(xi^2)=%.2f sum(xi*yi)=%.2f\n",
81 sum_xi, sum_yi, sum_xi_2, sum_xi_yi); */
83 *m = ( (double)sum_n * sum_xi_yi - sum_xi * sum_yi ) /
84 ( (double)sum_n * sum_xi_2 - sum_xi * sum_xi );
85 *b = (sum_yi / (double)sum_n) - (*m) * (sum_xi / (double)sum_n);
87 /* printf("slope = %.2f intercept = %.2f\n", *m, *b); */
92 return the least squares error:
98 double least_squares_error(double *x, double *y, int n, double m, double b) {
104 for ( i = 0; i < n; i++ ) {
105 error = y[i] - (m * x[i] + b);
106 sum += error * error;
107 // printf("%.2f %.2f\n", error, sum);
110 return ( sum / (double)n );
115 return the maximum least squares error:
119 double least_squares_max_error(double *x, double *y, int n, double m, double b){
121 double error, max_error;
125 for ( i = 0; i < n; i++ ) {
126 error = y[i] - (m * x[i] + b);
127 error = error * error;
128 if ( error > max_error ) {
133 return ( max_error );