least aquare method

hi,i would like to ask on how to fit curve with equation y=-ax-bx^2 using least square method.
The stochastic equation for (xi, yi) becomes:
1
2
3
yi = - a * xi – b * xi^2 + ei; 
//OR
ei = yi + a * xi + b * xi^2; 


where ei is the residual. Now:
(a) sum the squared residual over all (xi, yi)
(b) take the partial differentiations of this sum w.r.t a and b, which should be = 0 as per first-order conditions
(c) some assumptions about the distribution of the error (e.g. expectation 0 and orthogonal to x, y etc) will simplify the above equations further yielding the 'normal equations' which you can now solve for a^ and b^ which are estimates of a and b in the population
though your equation is not the classical OLS of econometrics any half-decent text-book on the subject (try the one by W H Greene if you will) should get you going
finally where C++ can come in handy is to arrange the X, Y data into containers like std::vector and then manipulating these vectors in solving the normal equations
you could also try Maximum Likelihood estimation using the <random> library if you wish to consider an alternative approach
Set:
squared error, S = Σ(y+ax+bx2)2

Set the partial derivatives to zero:
∂S/∂a = 2Σx(y+ax+bx2) = 0,
∂S/∂b = 2Σx2(y+ax+bx2) = 0
and solve the resulting simultaneous equations for a and b:
aΣx2+bΣx3=-Σxy
aΣx3+bΣx4=-Σx2y


In the code below you will want the sums on lines 23-31 and the expressions for a and b on lines 34-36.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
#include <iostream>
#include <iomanip>
#include <vector>
#include <fstream>
using namespace std;

int main()
{
   vector<double> x, y, yfit;

   // Read data
   double xval, yval;
   ifstream in( "data.txt" );
   while( in >> xval >> yval )
   {
      x.push_back( xval );
      y.push_back( yval );
   }
   in.close();

   // Compute sums
   int n = x.size();
   double s2 = 0, s3 = 0, s4 = 0, s1y = 0, s2y = 0;
   for ( int i = 0; i < n; i++ )
   {
      s2 += x[i] * x[i];
      s3 += x[i] * x[i] * x[i];
      s4 += x[i] * x[i] * x[i] * x[i];
      s1y += x[i] * y[i];
      s2y += x[i] * x[i] * y[i];
   }

   // Calculate fitted coefficients
   double denominator = s4 * s2 - s3 * s3;
   double a = ( s3 * s2y - s4 * s1y ) / denominator;
   double b = ( s3 * s1y - s2 * s2y ) / denominator;
   cout << "a = " << a << '\n';
   cout << "b = " << b << '\n';
   cout << '\n';

   // Fit data
   for ( int i = 0; i < n; i++ ) yfit.push_back( -a * x[i] - b * x[i] * x[i] );
 
   // Output tabulated data and curve fit
   #define SP << setw( 12 ) << setprecision( 4 ) << fixed <<
   cout SP "x" SP "y" SP "yfit" << '\n';
   for ( int i = 0; i < n; i++ ) cout SP x[i] SP y[i] SP yfit[i] << '\n';
}



data.txt
-2 -8
-1 -1
 0  0
 1 -5
 2 -16



Output:
a = 2
b = 3

           x           y        yfit
     -2.0000     -8.0000     -8.0000
     -1.0000     -1.0000     -1.0000
      0.0000      0.0000     -0.0000
      1.0000     -5.0000     -5.0000
      2.0000    -16.0000    -16.0000
Topic archived. No new replies allowed.