DEM/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 program is free software; you can redistribute it and/or modify
   8 // it under the terms of the GNU General Public License as published by
   9 // the Free Software Foundation; either version 2 of the License, or
  10 // (at your option) any later version.
  11 //
  12 // This program 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
  15 // GNU 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., 675 Mass Ave, Cambridge, MA 02139, USA.
  20 //
  21 // $Id$
  22 // (Log is kept at end of this file)
  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 void least_squares(double *x, double *y, int n, double *m, double *b) {
  45     double sum_xi, sum_yi, sum_xi_2, sum_xi_yi;
  46     int i;
  47
  48     sum_xi = sum_yi = sum_xi_2 = sum_xi_yi = 0.0;
  49
  50     for ( i = 0; i < n; i++ ) {
  51         sum_xi += x[i];
  52         sum_yi += y[i];
  53         sum_xi_2 += x[i] * x[i];
  54         sum_xi_yi += x[i] * y[i];
  55     }
  56
  57     /* printf("sum(xi)=%.2f  sum(yi)=%.2f  sum(xi^2)=%.2f  sum(xi*yi)=%.2f\n",
  58            sum_xi, sum_yi, sum_xi_2, sum_xi_yi); */
  59
  60     *m = ( (double)n * sum_xi_yi - sum_xi * sum_yi ) /
  61         ( (double)n * sum_xi_2 - sum_xi * sum_xi );
  62     *b = (sum_yi / (double)n) - (*m) * (sum_xi / (double)n);
  63
  64     /* printf("slope = %.2f  intercept = %.2f\n", *m, *b); */
  65 }
  66
  67
  68 /*
  69   return the least squares error:
  70
  71               (y[i] - y_hat[i])^2
  72               -------------------
  73                       n
  74  */
  75 double least_squares_error(double *x, double *y, int n, double m, double b) {
  76     int i;
  77     double error, sum;
  78
  79     sum = 0.0;
  80
  81     for ( i = 0; i < n; i++ ) {
  82         error = y[i] - (m * x[i] + b);
  83         sum += error * error;
  84         // printf("%.2f %.2f\n", error, sum);
  85     }
  86
  87     return ( sum / (double)n );
  88 }
  89
  90
  91 /*
  92   return the maximum least squares error:
  93
  94               (y[i] - y_hat[i])^2
  95  */
  96 double least_squares_max_error(double *x, double *y, int n, double m, double b){
  97     int i;
  98     double error, max_error;
  99
 100     max_error = 0.0;
 101
 102     for ( i = 0; i < n; i++ ) {
 103         error = y[i] - (m * x[i] + b);
 104         error = error * error;
 105         if ( error > max_error ) {
 106             max_error = error;
 107         }
 108     }
 109
 110     return ( max_error );
 111 }
 112
 113
 114 // $Log$
 115 // Revision 1.2  1998/04/21 17:03:41  curt
 116 // Prepairing for C++ integration.
 117 //
 118 // Revision 1.1  1998/04/08 22:57:24  curt
 119 // Adopted Gnu automake/autoconf system.
 120 //
 121 // Revision 1.1  1998/03/19 02:54:47  curt
 122 // Reorganized into a class lib called fgDEM.
 123 //
 124 // Revision 1.1  1997/10/13 17:02:35  curt
 125 // Initial revision.
 126 //