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.1  1998/04/08 22:57:24  curt
 116 /* Adopted Gnu automake/autoconf system.
 117 /*
 118  * Revision 1.1  1998/03/19 02:54:47  curt
 119  * Reorganized into a class lib called fgDEM.
 120  *
 121  * Revision 1.1  1997/10/13 17:02:35  curt
 122  * Initial revision.
 123  *
 124  */