src/least_squares_fitting/least_square.c

   1 //
   2 //     to fit a polynomial through data
   3 //                                           Michel Vallieres
   4 //                                           Spring 2002
   5 //
   6
   7
   8 #include <stdio.h>
   9 #include <math.h>
  10 #include <strings.h>
  11 #include <stdlib.h>
  12
  13
  14 //
  15 // standard deviation
  16 //
  17 int std_dev( double *residuals, double *sigma, int N )
  18 {
  19   int k;
  20   double sum;
  21
  22   sum = 0.0;
  23   for ( k=0 ; k<N ; k++ )
  24     {
  25       sum = sum + residuals[k]*residuals[k];
  26     }
  27   *sigma = sqrt( sum / (double)(N-1) );
  28 }
  29
  30
  31 //
  32 // evaluate fit function
  33 //
  34 int fit_function(
  35             double *alphas, double *x, double *y, double *fit,
  36             double *residuals, int N, char *fname )
  37 {
  38   int k;
  39   FILE *fp;
  40                               // function evaluation
  41   for ( k=0 ; k<N ; k++ )
  42     {
  43       fit[k] = alphas[0] + alphas[1]*x[k] + alphas[2]*x[k]*x[k];
  44       residuals[k] = fit[k]-y[k];
  45     }
  46                               // write fit function to file
  47   if ( ( fp = fopen( fname, "w" ) )  == NULL )
  48     {
  49        printf( " error in opening file (data_fitted) \n" );
  50        exit(1);
  51     }
  52   for ( k=0; k<N ; k++ )
  53     fprintf( fp, "%f %f %f %f\n",
  54                     x[k], y[k], fit[k], residuals[k] );
  55   close(fp);
  56 }
  57
  58
  59 //
  60 // alpha coefficients for the linear fit
  61 //
  62 int calculate_alphas_linear( double *alphas, double sums[5][5], int N )
  63 {
  64   double den;
  65
  66   den = N*sums[2][0] - sums[1][0]*sums[1][0];
  67   alphas[0] = ( sums[0][1]*sums[2][0] - sums[1][0]*sums[1][1] ) / den;
  68   alphas[1] = ( N*sums[1][1] - sums[1][0]*sums[0][1] ) / den;
  69   alphas[2] = 0.0;         // fit is linear
  70 }
  71
  72
  73 //
  74 // alpha coefficients for the polynomial fit
  75 //
  76 int calculate_alphas_square( double *alphas, double sums[5][5], int N )
  77 {
  78   double den;
  79
  80   den = - N*sums[4][0]*sums[2][0] + N*sums[3][0]*sums[3][0]
  81         + pow(sums[2][0],3) -  2.0*sums[3][0]*sums[2][0]*sums[1][0]
  82         + sums[4][0]*sums[1][0]*sums[1][0];
  83
  84   alphas[0] =   - sums[0][1]*sums[4][0]*sums[2][0]
  85                 + sums[0][1]*sums[3][0]*sums[3][0]
  86                 + sums[1][0]*sums[4][0]*sums[1][1]
  87                 - sums[1][0]*sums[3][0]*sums[2][1]
  88                 - sums[2][0]*sums[3][0]*sums[1][1]
  89                 + sums[2][0]*sums[2][0]*sums[2][1];
  90
  91   alphas[1] =   - N*sums[4][0]*sums[1][1]
  92                 + N*sums[3][0]*sums[2][1]
  93                 + sums[2][0]*sums[2][0]*sums[1][1]
  94                 - sums[2][0]*sums[1][0]*sums[2][1]
  95                 + sums[4][0]*sums[1][0]*sums[0][1]
  96                 - sums[3][0]*sums[2][0]*sums[0][1];
  97
  98   alphas[2] =   + N*sums[3][0]*sums[1][1]
  99                 - N*sums[2][0]*sums[2][1]
 100                 - sums[3][0]*sums[1][0]*sums[0][1]
 101                 - sums[2][0]*sums[1][0]*sums[1][1]
 102                 + sums[1][0]*sums[1][0]*sums[2][1]
 103                 + sums[2][0]*sums[2][0]*sums[0][1];
 104
 105   printf( "\n");
 106   printf( "  Denominator (check)  %f \n", den );
 107   printf( "  Numerator (alpha[0]) %f \n", alphas[0] );
 108   printf( "  Numerator (alpha[1]) %f \n", alphas[1] );
 109   printf( "  Numerator (alpha[2]) %f \n\n", alphas[2] );
 110
 111   alphas[0] = alphas[0] / den;
 112   alphas[1] = alphas[1] / den;
 113   alphas[2] = alphas[2] / den;
 114 }
 115
 116
 117 //
 118 // calculates all sums  needed in the fit formulae
 119 //
 120 int all_sums( double *x, double *y, double sums[5][5], int N )
 121 {
 122   int i,j;
 123   int k;
 124
 125   for ( i=0 ; i<5 ; i++ )
 126     for ( j=0 ; j<5 ; j++ )
 127       sums[i][j] = 0.0;
 128
 129   for ( i=0 ; i<5 ; i++ )
 130     for ( j=0 ; j<5 ; j++ )
 131       for ( k=0 ; k<N ; k++ )
 132         sums[i][j] = sums[i][j] + pow( x[k], i ) * pow( y[k], j );
 133 }
 134
 135
 136 //
 137 // main program to control the overall calculation
 138 //
 139 int main( int argc, char *argv[] )
 140 {
 141   int    i, j;
 142   int    N, k;
 143   FILE   *fp;
 144   double *x, *y, *fit, *residuals;
 145   double sums[5][5], alphas[5];
 146   char   fname[100];
 147   double sigma;
 148
 149                               // open data file
 150   if ( ( fp = fopen( "data_to_keep", "r" ) )  == NULL )
 151     {
 152        printf( " error in opening file (data_to_keeep) \n" );
 153        exit(1);
 154     }
 155                               // read in N ( # of data points )
 156   fscanf( fp, "%d\n", &N );
 157   printf( " Number of numbers: %d \n", N );
 158
 159                               // make space for arrays
 160   x = malloc( sizeof(double) * N );
 161   y = malloc( sizeof(double) * N );
 162   fit = malloc( sizeof(double) * N );
 163   residuals = malloc( sizeof(double) * N );
 164
 165                               // read in data
 166   for ( k=0; k<N ; k++ )
 167     fscanf( fp, "%lf %lf\n", &x[k], &y[k] );
 168
 169                               // close file
 170   fclose( fp );
 171
 172                               // perform all sums
 173   all_sums( x, y, sums, N );
 174                               // print sums matrix (check)
 175   printf( "Sums matrix\n" );
 176   for ( i=0 ; i<5 ; i++ )
 177     {
 178       for ( j=0 ; j<5 ; j++ )
 179         printf( " %f", sums[i][j] );
 180       printf(  "\n" );
 181     }
 182   printf(  "\n" );
 183
 184                               // Linear Fit
 185                               // ----------
 186
 187                               // calculate alphas
 188   calculate_alphas_linear( alphas, sums, N );
 189   printf( " Alphas: %f %f \n",
 190                          alphas[0], alphas[1] );
 191
 192                               // calculate the fit function
 193   strcpy( fname, "data_fitted_1" );
 194   fit_function( alphas, x, y, fit, residuals, N, fname );
 195
 196                               // standard deviation
 197   std_dev( residuals, &sigma, N );
 198   printf(" Standard deviation: %f\n", sigma );
 199
 200
 201                               // Quadratic Fit
 202                               // -------------
 203
 204                               // calculate alphas
 205   calculate_alphas_square( alphas, sums, N );
 206   printf( " Alphas: %f %f %f\n",
 207                          alphas[0], alphas[1], alphas[2] );
 208
 209                               // calculate the fit function
 210   strcpy( fname, "data_fitted_2" );
 211   fit_function( alphas, x, y, fit, residuals, N, fname );
 212
 213                               // standard deviation
 214   std_dev( residuals, &sigma, N );
 215   printf(" Standard deviation: %f\n", sigma );
 216
 217
 218                               // free memory
 219   free( x );
 220   free( y );
 221   free( fit );
 222                               // done
 223   exit(0);
 224 }
 225
 226
 227
 228
 229
 230
 231
 232
 233