1 // leastsqs.c -- Implements a simple linear least squares best fit routine
3 // Written by Curtis Olson, started September 1997.
5 // Copyright (C) 1997 Curtis L. Olson - curt@infoplane.com
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.
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.
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.
22 // (Log is kept at end of this file)
28 #include "leastsqs.hxx"
36 n*sum(xi*yi) - (sum(xi)*sum(yi))
37 b1 = --------------------------------
38 n*sum(xi^2) - (sum(xi))^2
41 b0 = sum(yi)/n - b1*(sum(xi)/n)
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;
48 sum_xi = sum_yi = sum_xi_2 = sum_xi_yi = 0.0;
50 for ( i = 0; i < n; i++ ) {
53 sum_xi_2 += x[i] * x[i];
54 sum_xi_yi += x[i] * y[i];
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); */
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);
64 /* printf("slope = %.2f intercept = %.2f\n", *m, *b); */
69 return the least squares error:
75 double least_squares_error(double *x, double *y, int n, double m, double b) {
81 for ( i = 0; i < n; i++ ) {
82 error = y[i] - (m * x[i] + b);
84 // printf("%.2f %.2f\n", error, sum);
87 return ( sum / (double)n );
92 return the maximum least squares error:
96 double least_squares_max_error(double *x, double *y, int n, double m, double b){
98 double error, max_error;
102 for ( i = 0; i < n; i++ ) {
103 error = y[i] - (m * x[i] + b);
104 error = error * error;
105 if ( error > max_error ) {
110 return ( max_error );
115 // Revision 1.2 1998/04/21 17:03:41 curt
116 // Prepairing for C++ integration.
118 // Revision 1.1 1998/04/08 22:57:24 curt
119 // Adopted Gnu automake/autoconf system.
121 // Revision 1.1 1998/03/19 02:54:47 curt
122 // Reorganized into a class lib called fgDEM.
124 // Revision 1.1 1997/10/13 17:02:35 curt