From 7b22b8cd929d1dded5b678076875d393c27e9d3d Mon Sep 17 00:00:00 2001 From: curt Date: Sat, 13 Mar 1999 17:34:44 +0000 Subject: [PATCH] Moved to math subdirectory. --- Math/Makefile.am | 1 + Math/leastsqs.cxx | 152 ++++++++++++++++++++++++++++++++++++++++++++++ Math/leastsqs.hxx | 90 +++++++++++++++++++++++++++ 3 files changed, 243 insertions(+) create mode 100644 Math/leastsqs.cxx create mode 100644 Math/leastsqs.hxx diff --git a/Math/Makefile.am b/Math/Makefile.am index e6654827..eda2bc00 100644 --- a/Math/Makefile.am +++ b/Math/Makefile.am @@ -8,6 +8,7 @@ libMath_a_SOURCES = \ fg_geodesy.cxx fg_geodesy.hxx \ fg_random.c fg_random.h \ interpolater.cxx interpolater.hxx \ + leastsqs.cxx leastsqs.hxx \ mat3.h mat3defs.h mat3err.h \ point3d.hxx \ polar3d.cxx polar3d.hxx \ diff --git a/Math/leastsqs.cxx b/Math/leastsqs.cxx new file mode 100644 index 00000000..28e2bd35 --- /dev/null +++ b/Math/leastsqs.cxx @@ -0,0 +1,152 @@ +// leastsqs.c -- Implements a simple linear least squares best fit routine +// +// Written by Curtis Olson, started September 1997. +// +// Copyright (C) 1997 Curtis L. Olson - curt@infoplane.com +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. +// +// $Id$ +// (Log is kept at end of this file) +// + + +#include + +#include "leastsqs.hxx" + + +/* +Least squares fit: + +y = b0 + b1x + + n*sum(xi*yi) - (sum(xi)*sum(yi)) +b1 = -------------------------------- + n*sum(xi^2) - (sum(xi))^2 + + +b0 = sum(yi)/n - b1*(sum(xi)/n) +*/ + +double sum_xi, sum_yi, sum_xi_2, sum_xi_yi; +int sum_n; + +void least_squares(double *x, double *y, int n, double *m, double *b) { + int i; + + sum_xi = sum_yi = sum_xi_2 = sum_xi_yi = 0.0; + sum_n = n; + + for ( i = 0; i < n; i++ ) { + sum_xi += x[i]; + sum_yi += y[i]; + sum_xi_2 += x[i] * x[i]; + sum_xi_yi += x[i] * y[i]; + } + + /* printf("sum(xi)=%.2f sum(yi)=%.2f sum(xi^2)=%.2f sum(xi*yi)=%.2f\n", + sum_xi, sum_yi, sum_xi_2, sum_xi_yi); */ + + *m = ( (double)sum_n * sum_xi_yi - sum_xi * sum_yi ) / + ( (double)sum_n * sum_xi_2 - sum_xi * sum_xi ); + *b = (sum_yi / (double)sum_n) - (*m) * (sum_xi / (double)sum_n); + + /* printf("slope = %.2f intercept = %.2f\n", *m, *b); */ +} + + +/* incrimentally update existing values with a new data point */ +void least_squares_update(double x, double y, double *m, double *b) { + ++sum_n; + + sum_xi += x; + sum_yi += y; + sum_xi_2 += x * x; + sum_xi_yi += x * y; + + /* printf("sum(xi)=%.2f sum(yi)=%.2f sum(xi^2)=%.2f sum(xi*yi)=%.2f\n", + sum_xi, sum_yi, sum_xi_2, sum_xi_yi); */ + + *m = ( (double)sum_n * sum_xi_yi - sum_xi * sum_yi ) / + ( (double)sum_n * sum_xi_2 - sum_xi * sum_xi ); + *b = (sum_yi / (double)sum_n) - (*m) * (sum_xi / (double)sum_n); + + /* printf("slope = %.2f intercept = %.2f\n", *m, *b); */ +} + + +/* + return the least squares error: + + (y[i] - y_hat[i])^2 + ------------------- + n + */ +double least_squares_error(double *x, double *y, int n, double m, double b) { + int i; + double error, sum; + + sum = 0.0; + + for ( i = 0; i < n; i++ ) { + error = y[i] - (m * x[i] + b); + sum += error * error; + // printf("%.2f %.2f\n", error, sum); + } + + return ( sum / (double)n ); +} + + +/* + return the maximum least squares error: + + (y[i] - y_hat[i])^2 + */ +double least_squares_max_error(double *x, double *y, int n, double m, double b){ + int i; + double error, max_error; + + max_error = 0.0; + + for ( i = 0; i < n; i++ ) { + error = y[i] - (m * x[i] + b); + error = error * error; + if ( error > max_error ) { + max_error = error; + } + } + + return ( max_error ); +} + + +// $Log$ +// Revision 1.1 1999/03/13 17:34:45 curt +// Moved to math subdirectory. +// +// Revision 1.2 1998/04/21 17:03:41 curt +// Prepairing for C++ integration. +// +// Revision 1.1 1998/04/08 22:57:24 curt +// Adopted Gnu automake/autoconf system. +// +// Revision 1.1 1998/03/19 02:54:47 curt +// Reorganized into a class lib called fgDEM. +// +// Revision 1.1 1997/10/13 17:02:35 curt +// Initial revision. +// diff --git a/Math/leastsqs.hxx b/Math/leastsqs.hxx new file mode 100644 index 00000000..d8b40c80 --- /dev/null +++ b/Math/leastsqs.hxx @@ -0,0 +1,90 @@ +// leastsqs.h -- Implements a simple linear least squares best fit routine +// +// Written by Curtis Olson, started September 1997. +// +// Copyright (C) 1997 Curtis L. Olson - curt@infoplane.com +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. +// +// $Id$ +// (Log is kept at end of this file) +/// + + +#ifndef _LEASTSQS_H +#define _LEASTSQS_H + + +#ifndef __cplusplus +# error This library requires C++ +#endif + + +/* +Least squares fit: + +y = b0 + b1x + + n*sum(xi*yi) - (sum(xi)*sum(yi)) +b1 = -------------------------------- + n*sum(xi^2) - (sum(xi))^2 + + +b0 = sum(yi)/n - b1*(sum(xi)/n) +*/ + +void least_squares(double *x, double *y, int n, double *m, double *b); + +/* incrimentally update existing values with a new data point */ +void least_squares_update(double x, double y, double *m, double *b); + + +/* + return the least squares error: + + (y[i] - y_hat[i])^2 + ------------------- + n +*/ +double least_squares_error(double *x, double *y, int n, double m, double b); + + +/* + return the maximum least squares error: + + (y[i] - y_hat[i])^2 +*/ +double least_squares_max_error(double *x, double *y, int n, double m, double b); + + +#endif // _LEASTSQS_H + + +// $Log$ +// Revision 1.1 1999/03/13 17:34:45 curt +// Moved to math subdirectory. +// +// Revision 1.2 1998/04/21 17:03:42 curt +// Prepairing for C++ integration. +// +// Revision 1.1 1998/04/08 22:57:25 curt +// Adopted Gnu automake/autoconf system. +// +// Revision 1.1 1998/03/19 02:54:48 curt +// Reorganized into a class lib called fgDEM. +// +// Revision 1.1 1997/10/13 17:02:35 curt +// Initial revision. +// -- 2.39.5