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