simgear/math/leastsqs.cxx

   1 // leastsqs.c -- Implements a simple linear least squares best fit routine
   2 //
   3 // Written by Curtis Olson, started September 1997.
   4 //
   5 // Copyright (C) 1997  Curtis L. Olson  - http://www.flightgear.org/~curt
   6 //
   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.
  11 //
  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.
  16 //
  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.
  20 //
  21 // $Id$
  22 //
  23
  24
  25 #include <stdio.h>
  26
  27 #include "leastsqs.hxx"
  28
  29
  30 /*
  31 Least squares fit:
  32
  33 y = b0 + b1x
  34
  35      n*sum(xi*yi) - (sum(xi)*sum(yi))
  36 b1 = --------------------------------
  37      n*sum(xi^2) - (sum(xi))^2
  38
  39
  40 b0 = sum(yi)/n - b1*(sum(xi)/n)
  41 */
  42
  43 double sum_xi, sum_yi, sum_xi_2, sum_xi_yi;
  44 int sum_n;
  45
  46 void least_squares(double *x, double *y, int n, double *m, double *b) {
  47     int i;
  48
  49     sum_xi = sum_yi = sum_xi_2 = sum_xi_yi = 0.0;
  50     sum_n = n;
  51
  52     for ( i = 0; i < n; i++ ) {
  53         sum_xi += x[i];
  54         sum_yi += y[i];
  55         sum_xi_2 += x[i] * x[i];
  56         sum_xi_yi += x[i] * y[i];
  57     }
  58
  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); */
  61
  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);
  65
  66     /* printf("slope = %.2f  intercept = %.2f\n", *m, *b); */
  67 }
  68
  69
  70 /* incrimentally update existing values with a new data point */
  71 void least_squares_update(double x, double y, double *m, double *b) {
  72     ++sum_n;
  73
  74     sum_xi += x;
  75     sum_yi += y;
  76     sum_xi_2 += x * x;
  77     sum_xi_yi += x * y;
  78
  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); */
  81
  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);
  85
  86     /* printf("slope = %.2f  intercept = %.2f\n", *m, *b); */
  87 }
  88
  89
  90 /*
  91   return the least squares error:
  92
  93               (y[i] - y_hat[i])^2
  94               -------------------
  95                       n
  96  */
  97 double least_squares_error(double *x, double *y, int n, double m, double b) {
  98     int i;
  99     double error, sum;
 100
 101     sum = 0.0;
 102
 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);
 107     }
 108
 109     return ( sum / (double)n );
 110 }
 111
 112
 113 /*
 114   return the maximum least squares error:
 115
 116               (y[i] - y_hat[i])^2
 117  */
 118 double least_squares_max_error(double *x, double *y, int n, double m, double b){
 119     int i;
 120     double error, max_error;
 121
 122     max_error = 0.0;
 123
 124     for ( i = 0; i < n; i++ ) {
 125         error = y[i] - (m * x[i] + b);
 126         error = error * error;
 127         if ( error > max_error ) {
 128             max_error = error;
 129         }
 130     }
 131
 132     return ( max_error );
 133 }
 134
 135