Main Page   Namespace List   Class Hierarchy   Compound List   File List   Compound Members   File Members  

lregress.cpp

Go to the documentation of this file.
00001 
00002 #include <math.h>
00003 
00004 #include "lregress.h"
00005 
00006 using namespace std;
00007 
00013 double lregress::meanX( vector<lregress::point> &data ) const
00014 {
00015   double mean = 0.0;
00016   double sum = 0.0;
00017 
00018   const size_t N = data.size();
00019   size_t i;
00020   for (i = 0; i < N; i++) {
00021     sum = sum + data[i].x();
00022   }
00023   mean = sum / N;
00024   return mean;
00025 } // meanX
00026 
00032 double lregress::meanY( vector<lregress::point> &data ) const
00033 {
00034   double mean = 0.0;
00035   double sum = 0.0;
00036 
00037   const size_t N = data.size();
00038   size_t i;
00039   for (i = 0; i < N; i++) {
00040     sum = sum + data[i].y();
00041   }
00042   mean = sum / N;
00043   return mean;
00044 } // meanY
00045 
00046 
00047 
00064 void lregress::lr( vector<lregress::point> &data,
00065                    lregress::lineInfo &info ) const
00066 {
00067   const size_t N = data.size();
00068   if (N > 0) {
00069 
00070     double muX = meanX( data );
00071 
00072     double muY = meanY( data );
00073 
00074     //     N-1
00075     //     ---
00076     //     \   (Xi - meanX)(Yi - meanY)
00077     //     /__ 
00078     //     i=0
00079     // b = -----------------------------
00080     //     N-1
00081     //     ---             2
00082     //     \   (Xi - meanX)
00083     //     /__ 
00084     //     i=0
00085     //
00086 
00087     double SSxy = 0;
00088     double SSxx = 0;
00089     double SSyy = 0;
00090     double Sx = 0;
00091     double Sy = 0;
00092     double Sxy = 0;
00093     double SSy = 0;
00094     double SSx = 0;
00095     for (size_t i = 0; i < N; i++) {
00096       Sx = Sx + data[i].x();
00097       Sy = Sy + data[i].y();
00098       Sxy = Sxy + (data[i].x() * data[i].y());
00099       SSx = SSx + (data[i].x() * data[i].x());
00100       SSy = SSy + (data[i].y() * data[i].y());
00101       double subX = (data[i].x() - muX);
00102       double subY = (data[i].y() - muY);
00103       SSyy = SSyy + subY * subY;
00104       SSxy = SSxy + subX * subY;
00105       SSxx = SSxx + subX * subX;
00106     }
00107     
00108     // slope
00109     double b = SSxy / SSxx;
00110     // intercept
00111     double a = muY - b * muX;
00112 
00113     // standard deviation of the points
00114     double stddevPoints = sqrt( (SSyy - b * SSxy)/(N-2) );
00115 
00116     // Error of the slope
00117     double bError = stddevPoints / sqrt( SSxx );
00118 
00119     double r2Numerator = (N * Sxy) - (Sx * Sy);
00120     double r2Denominator = ((N*SSx) - (Sx * Sx))*((N*SSy) - (Sy * Sy));
00121     double r2 = (r2Numerator * r2Numerator) / r2Denominator;
00122 
00123     double signB = (b < 0) ? -1.0 : 1.0;
00124 
00125     double r = signB * sqrt( r2 );
00126   
00127     info.corr( r );
00128     info.stddev( stddevPoints );
00129     info.slopeErr( bError );
00130     info.slope( b );
00131     info.intercept( a );
00132   } // if N > 0
00133 
00134 } // lineInfo

Generated at Tue May 27 21:56:16 2003 for Wavelet compression, determinism and time series forecasting by doxygen1.2.8.1 written by Dimitri van Heesch, © 1997-2001