Best Line(Curve) Fitting - Least Square Method

          Curve fitting is very fun with our own coding that fits a best curve for a given set of points. Although it is fun, it is also an essential friend for students for plotting their readings. Most of the engineering stream students have to submit their observation with well plotted graph. But what they actually do is draw the graph and adjust the readings according to the their graph. I hope this program will help them to plot their graph more appropriately.


Theory
        The name least square denote the squares of the errors are the least possible values. What are the errors? The difference between our actual reading and the value of the chosen equation.
        Let us assume we are going to fit a straight line for the given set of readings. The equation of the line is : y = mx+c. This equation has two parameters x and y thus it has it own (x,y) values and we have our own set of (x,y) values observed from the experiment. The difference between these y (from equation) and y (from experiment) is called error.


Our Expt: Ohm's Law.

The circuit on the left is to prove the Ohm's Law V = IR.
Definition: Ohm's law states that the current through a conductor between two points is directly proportional to the potential difference across the two points, and inversely proportional to the resistance between them.


Assume that this is the set of readings we have:    

X ( Current )	Y (Voltage )
1	2.3
2	4.6
3	5.8
4	7.5
5	9.8

                                            

        This plot is not a perfect straight line. So we try to fit a straight line such that it is nearer to all of the points. What we are going to do is place a scale on the graph sheet and adjust it until it becomes nearer to all points and then we draw the line. Mathematically we fit a line y = mx+c, then we find the values of 'm' and 'c' so that it resemble the line we draw using the experimental values. For more theory google for it.

Coding:          

       We use the matrix library that was discussed in my previous post

File: lsm.h
   #include<matrix.h>
   typedef struct readings {
     float *x;
     float *y;
     int len;
   }Readings;
   typedef struct lsm_line {
     float *x;
     float *y;
     float *x2;
     float *xy;
     float X , Y , XY , X2;
     float *newy;
     float Dev;
     int len;
   }LsmLine;
   Readings* lsm_new_Readings(int n);
   LsmLine* lsm_new_LsmLine(Readings *Read);
   int lsm_del_Readings(Readings *anchor);
   int lsm_Read_farray(float *array,int len);
   int lsm_Print_farray(float *array,int len);
   int lsm_Calc_LsmLine(LsmLine* line);
   int lsm_Print_LsmLine(LsmLine* line); 
File : lsm.c

     #include "lsm.h"

     Readings* lsm_new_Readings(int len)

     {

        Readings *newRead = (Readings*)malloc(sizeof(Readings));

        newRead->x = (float*)malloc(sizeof(float)*len);

        newRead->y = (float*)malloc(sizeof(float)*len);

        newRead->len = len;

        return newRead;

     }

     LsmLine* lsm_new_LsmLine(Readings *read)

     {

        LsmLine *newLine = (LsmLine*)malloc(sizeof(LsmLine));

        int len = read->len;

        newLine->x = (float*)malloc(sizeof(float)*len);

        newLine->y = (float*)malloc(sizeof(float)*len);

        newLine->xy = (float*)malloc(sizeof(float)*len);

        newLine->x2= (float*)malloc(sizeof(float)*len);

        newLine->newy= (float*)malloc(sizeof(float)*len);

        newLine->len = len;

        newLine->X  = 0;

        newLine->Y  = 0;

        newLine->XY = 0;

        newLine->X2 = 0;

        newLine->Dev=0;

        int i;

        newLine->len = read->len;

        for(i = 0;i<newLine->len;i++){

            newLine->x[i]=read->x[i];

            newLine->y[i]=read->y[i];

            newLine->xy[i]=read->x[i]*read->y[i];

            newLine->x2[i]=read->x[i]*read->x[i];

            newLine->X +=  newLine->x[i];

            newLine->Y +=  newLine->y[i];

            newLine->XY +=  newLine->xy[i];

            newLine->X2 +=  newLine->x2[i];

        }

        return newLine;

  }

  
  int lsm_Calc_LsmLine(LsmLine *Line)

  {

     int i;

     Matrix *lhs = matrix_new_Matrix(2,2);

     Matrix *rhs = matrix_new_Matrix(2,1);

     Assuming Eqn is y=ax+b

     lhs->mat[0][0] = Line->X;

     lhs->mat[0][1] = Line->len;

     lhs->mat[1][0] = Line->X2;

     lhs->mat[1][2] = Line->X;

     rhs->mat[0][0]=Line->Y;

     rhs->mat[1][0]=Line->XY;

     MatrixEqn *Eqn = matrix_new_MatrixEqn(lhs,rhs);

     matrix_Solve_Det(Eqn);

     float a = Eqn->Soln->mat[0][0];

     float b = Eqn->Soln->mat[1][0] ;

     for(i=0;i<Line->len;i++){

       Line->newy[i] =  a*Line->x[i]+b;

       Line->Dev += (Line->y[i]-Line->newy[i])*(Line->y[i]-Line->newy[i]);

     }

 }

 int lsm_Print_LsmLine(LsmLine* line)

 {

    printf("\nS.No\tX\t\tY\t\tX2\t\tXY\t\tNewY\n");

    int i;

    for (i=0;i<line->len;i++){

       printf("\n%d",i+1);

       printf("\t%0.3f",line->x[i]);

       printf("\t\t%0.3f",line->y[i]);

       printf("\t\t%0.3f",line->x2[i]);

       printf("\t%0.3f",line->xy[i]);

       printf("\t\t%0.3f",line->newy[i]);

    }

    printf("\n\nTotal");

    printf("\t%0.3f",line->X);

    printf("\t\t%0.3f",line->Y);

    printf("\t\t%0.3f",line->X2);

    printf("\t%0.3f",line->XY);

    printf("\n\n%0.3f = %0.3fa + %db ",line->Y,line->X,line->len);

    printf("\n%0.3f = %0.3fa + %0.3fb ",line->XY,line->X2,line->X);

    printf("\n\nDeviation = %0.3f",line->Dev);

    return 0;
  }

  
  void lsm_Read_farray(float *array,int len)

  {

     int i=0;

     for(i=0;i<len;i++){

        scanf("%f",&array[i]);

     }

     return 0;

  }

  
   lsm_Print_farray(float *array,int len)

   {

      int i=0;

      for(i=0;i<len;i++){

         printf("%0.3f",array[i]);

      }

      return 0;

   }
Sample program to test.

       #include<stdio.h>

   #include"lsm.h"

   int main()

   {

       printf("\n\n\t\t\tLEAST SQUARE METHOD FOR BEST LINE");

       int n,choice;

       printf("\n\nEnter the number of readings\t:");

       scanf("%d",&n);

       Readings *read = lsm_new_Readings(n);

       printf("\nEnter the type of Variation of X:");

       printf("\n\t1.Constant Steps of X");

       printf("\n\t2.Variable Steps of X\n");

       printf("\nEnter the Choice\t:");

       scanf("%d",&choice);

       switch(choice){

         case 1:   printf("\nEnter First : Step : Last of X\t:") ;

                   float first , step , last;

                   scanf("%f %f %f",&first,&step,&last);

                   int i;

                   float value;

                   for( i=0,value=first;value<=last; value +=step,i++){

                       read->x[i] = value;

                   }

                   break;

         case 2:   printf("\nEnter X\t:\n");

                   lsm_Read_farray(read->x,n);

                   break;

         default: printf("\nWrong Option Restart the program!!!");

                   break;

      }

      printf("\nEnter Y\t:\n");

      lsm_Read_farray(read->y,n);

      LsmLine *line = lsm_new_LsmLine(read);

      lsm_Calc_LsmLine(line);

      lsm_Print_LsmLine(line);

      return 0;

}

         We can extend these project to various curves, and choose the curve which has least summation of error 'S', since all set of points will not be as above.

Matrix Manipulation Using C

            Matrix  is a great invention of Mathematics. It is very useful in solving system of equations through computer. Of course famous numerical computing platform, MATLAB is the well known application of matrix.
            In MATLAB everything is matrix. The images are also represented as a matrix in ease. It is easy to perform operations over matrices. This is why the matrix manipulation is an important.

Download the code here

            Here is my matrix library. It may be inefficient, since it is indented for educational purpose. I want to keep things simple so that the beginner also be able to understand. This article is also helpful for learning memory allocation in C, Pointer and its applications, recursion,  modularity.
                                                                
CodeListing

1. /*----------------------------------------*/
2. /* Program: matrix.h                      */
3. /* Author: Vanangamudi                    */
4. /* Date: 15/05/2011                       */
5. /*----------------------------------------*/
6. #include<stdio.h>
7. #ifndef _MATRIX_H_
8. #define _MATRIX_H_
9. 
   
10. typedef struct matrix {
11. float **mat;
12. int m,n;
13. }Matrix;
14. 
    
15. typedef struct matrixEquation {
16. Matrix *lhs;
17. Matrix *rhs;
18. Matrix *Soln;
19.  }MatrixEqn;
20. 
    
21. Matrix * matrix_new_Matrix(int m,int n);
22. MatrixEqn * matrix_new_MatrixEqn(Matrix* lhs,Matrix *rhs);
23. 
    
24. float matrix_Det(Matrix *matA);
25. Matrix* matrix_Cofactor(Matrix *matA);
26. Matrix* matrix_Minor(Matrix *matA, int m, int n);
27. Matrix* matrix_Adjoint(Matrix *matA);
28. Matrix* matrix_Inverse(Matrix *matA);
29. 
    
30. int matrix_Solve_Det(MatrixEqn *matEqn);
31. int matrix_Solve_Inv(MatrixEqn *matEqn);
32. 
    
33. Matrix* matrix_Mult(Matrix *matA , Matrix *matB);
34. Matrix* matrix_Add(Matrix *matA , Matrix *matB);
35. Matrix* matrix_Sub(Matrix *matA , Matrix *matB);
36. Matrix* matrix_tran(Matrix *matA);
37. Matrix* matrix_ScalarMultiply(float X , Matrix *matA);
38. Matrix* matrix_Duplicate(Matrix *matA);
39. Matrix* matrix_ColumnReplace(Matrix *matA, Matrix* matB ,int IndA,int IndB );
40. 
    
41. int matrix_Print(Matrix *matA);
42. int matrix_Read(Matrix *matA);
43. 
    
44. #endif       //_MATRIX_H_

1. /*-----------------------------------------------*/
2. /* Program: matrix.c                             */
3. /* Author: Vanangamudi                           */
4. /* Date: 15/05/2011                              */
5. /*-----------------------------------------------*/
6. 
   
7. #include "matrix.h"
8. #include<malloc.h>
9. #include<stdio.h>
10. 
    
11. //Constructors
12. Matrix* matrix_new_Matrix(int m, int n)
13. {
14.        Matrix *newMat = (Matrix*)malloc(sizeof(float));
15.        newMat->mat = (float*)malloc(sizeof(float)*m);
16.     int i = 0;
17.     for(i =0 ; i < m ; i++){
18.           newMat->mat[i] = (float*)malloc(sizeof(float)*n);
19.     }
20.       newMat->m = m;
21.       newMat->n = n;
22.     return newMat;
23. }
24. 
    
25. MatrixEqn * matrix_new_MatrixEqn(Matrix *lhs,Matrix *rhs)
26. {
27.       MatrixEqn* newEqn = (MatrixEqn*)malloc(sizeof(float));
28.       newEqn->lhs = matrix_Duplicate(lhs); 
29.       newEqn->rhs = matrix_Duplicate(rhs);
30.       newEqn->Soln = matrix_new_Matrix(rhs->m,1);
31.     return newEqn;
32. }
33. 
    
34. // Matrix Operations
35. Matrix* matrix_Minor(Matrix *matA, int m, int n)
36. {
37.     int i,j;
38.     int p,q;
39.     Matrix *minor = matrix_new_Matrix(matA->m-1,matA->n-1);
40.     for(i = 0,p = 0; i < matA->m , p < minor->m; i++){
41.       if( i != m-1){
42.          for(j=0, q=0; j <matA->n , q < minor->n; j++){
43.             if(j != n-1){
44.               minor->mat[p][q] = matA->mat[i][j];
45.               q++;
46.             }
47.          }
48.          p++;
49.       }
50.     }
51.     return minor;
52. }
53. 
    
54. float matrix_Det(Matrix *matA)
55. {
56.     float det=0;
57.     if((matA->m==1)&&(matA->n==1)){
58.           det = matA->mat[0][0];
59.     }else {
60.       int i;
61.       int sign=1;
62.       for(i = 0; i<matA->m ; i++){
63.          det += sign*matA->mat[0][i]*matrix_Det(matrix_Minor(matA,0,i));
64.           sign = -1*sign;
65.       }
66.     }
67.     return det;
68. }
69. 
    
70. Matrix* matrix_Cofactor(Matrix *matA)
71. {
72.     Matrix *Cofact = matrix_new_Matrix(matA->m,matA->n);
73.     int i,j;
74.     int sign =1;
75.     for(i =0 ; i<matA->m;i++){
76.        for(j=0 ; j<matA->n ; j++){
77.           Cofact->mat[i][j] = sign*matrix_Det(matrix_Minor(matA,i,j));
78.           sign = sign * (-1);
79.        }
80.        if(matA->m%2 == 0){
81.               sign = (-1)*sign;
82.        }
83.     }
84.     return Cofact;
85. }
86. 
    
87. Matrix* matrix_Adjoint(Matrix *matA)
88. {
89.      Matrix *Cofact = matrix_Cofactor(matA);
90.      Matrix *Adj = matrix_tran(Cofact);
91.      return Adj;
92. }
93. 
    
94. Matrix* matrix_Inverse(Matrix *matA)
95. {
96.     Matrix *Inv = matrix_Adjoint(matA);
97.     float Det = matrix_Det(matA);
98.     Inv = matrix_ScalarMultiply(1/Det, Inv);
99.     return Inv;
100. }
101. 
     
102. //Arithmetics
103. Matrix* matrix_Mult(Matrix *matA , Matrix *matB)
104. {
105.      if(matA->n != matB->m){
106.         printf("\nMatrices cannot be multiplied");
107.         return NULL;
108.      }else{
109.              Matrix *matC = matrix_new_Matrix(matA->m,matB->n);
110.         int i,j,k;
111.         for(i = 0 ;i< matB->m; i++){
112.            for(j=0; j < matB->n ; j++){
113.               matC->mat[i][j] = 0;
114.               for( k =0;k < matA->n;k++){
115.                  matC->mat[i][j]+=matA->mat[i][k]*matB->mat[k][j];
116.               }
117.             }
118.          }
119.         return matC;
120.      }
121. }
122. 
     
123. Matrix* matrix_ScalarMultiply(float X , Matrix *matA)
124. {
125.         Matrix *matB = matrix_new_Matrix(matA->m,matA->n);
126.      int i,j;
127.      for(i = 0 ;i< matA->m; i++){
128.         for(j=0; j < matA->n ; j++){
129.            matB->mat[i][j] = X * matA->mat[i][j];
130.         }
131.      }
132.      return matB;
133. }
134. 
     
135. Matrix* matrix_Add(Matrix *matA , Matrix *matB)
136. {
137.         Matrix *matC = matrix_new_Matrix(matA->m,matA->n);
138.      int i,j;
139.      for(i = 0 ;i< matA->m; i++){
140.         for(j=0; j < matA->n ; j++){
141.           matC->mat[i][j] = matA->mat[i][j]+matB->mat[i][j];
142.         }
143.      }
144.       return matC;
145. }
146. 
     
147. Matrix* matrix_Sub(Matrix *matA , Matrix *matB)
148. {
149.         Matrix *matC = matrix_new_Matrix(matA->m,matA->n);
150.      int i,j;
151.      for(i = 0 ;i< matA->m; i++){
152.         for(j=0; j < matA->n ; j++){
153.            matC->mat[i][j] = matA->mat[i][j]-matB->mat[i][j];
154.         }
155.      }
156.      return matC;
157. }
158. 
     
159. Matrix* matrix_tran(Matrix *matA)
160. {
161.     Matrix *matB = matrix_new_Matrix(matA->n,matA->m);
162.     int i,j;
163.     for(i = 0 ;i< matB->m; i++){
164.        for(j=0; j < matB->n ; j++){
165.           matB->mat[i][j] = matA->mat[j][i];
166.        }
167.     }
168.     return matB;
169. }
170. 
     
171. //Matrix Eqn Solvers
172. int matrix_Solve_Inv(MatrixEqn *matEqn)
173. {
174.     Matrix* Inv = matrix_Inverse(matEqn->lhs);
175.     matEqn->Soln = matrix_Mult(Inv,matEqn->rhs);
176. }
177. 
     
178. int matrix_Solve_Det(MatrixEqn *matEqn)
179. {
180.     float Det = matrix_Det(matEqn->lhs);
181.     Matrix* matA;
182.     float DetXYZ;
183.     int i;
184.     for(i = 0;i <matEqn->Soln->m;i++){
185.       matA = matrix_ColumnReplace(matEqn->lhs,matEqn->rhs,i+1,1);
186.       DetXYZ = matrix_Det(matA);
187.       matEqn->Soln->mat[i][0] = DetXYZ/Det;
188.     }
189. }
190. 
     
191. //Misc
192. Matrix* matrix_ColumnReplace(Matrix *matA, Matrix* matB ,int IndA,int IndB )
193. {
194.        Matrix* matC = matrix_Duplicate(matA);
195.     int i;
196.     for(i = 0;i<matA->n; i++){
197.             matC->mat[i][IndA-1] = matB->mat[i][IndB-1];
198.     }
199.     return matC;
200. }
201. 
     
202. Matrix* matrix_Duplicate(Matrix *matA)
203. {
204.        Matrix *matB = matrix_new_Matrix(matA->m,matA->n);
205.     int i,j;
206.     for(i = 0 ;i< matB->m; i++){
207.        for(j=0; j < matB->n ; j++){
208.                 matB->mat[i][j] = matA->mat[i][j];
209.        }
210.     }
211.     return matB;
212. }
213. 
     
214. //IO Operations
215. int matrix_Read(Matrix *matA)
216. {
217.      int i=0 , j=0;
218.         printf("\nEnter the Elements:\n");
219.      for(i = 0; i < matA->m; i++){
220.         for(j = 0; j < matA->n; j++ ){
221.                  scanf("%f",&matA->mat[i][j]);
222.         }
223.      }
224. }
225. 
     
226. int matrix_Print(Matrix *matA)
227. {
228.      int i=0 , j=0;
229.         printf("\n");
230.      for(i = 0; i < matA->m; i++){
231.         for(j = 0; j < matA->n; j++ ){
232.                  printf("%0.3f\t",matA->mat[i][j]);
233.         }
234.              printf("\n\n");
235.      }
236.      return 0;
237. }

Sample Program to Test

1. /*-----------------------------------------------------*/
2. /* Program: matrix_test.c                              */
3. /* Author  : Vanangamudi                               */
4. /* Date     :15/05/2011                                */
5. /*-----------------------------------------------------*/
6. 
   
7. #include "matrix.h"
8. #include<stdio.h>
9. 
   
10. int main()
11. {
12.     Matrix *matA  = matrix_new_Matrix(3,3);
13.     Matrix *matB  = matrix_new_Matrix(3,3);
14. 
    
15.     printf("\t\t|--------------------------|");
16.     printf("\n\t\t|**************************|");
17.     printf("\n\t\t|*** MATRIX ILLUSTRATION **|\n");
18.     printf("\t\t|**************************|\n");
19.     printf("\t\t|--------------------------|\n\n");
20. 
    
21.     matrix_Read(matA);
22.     matrix_Read(matB);
23. 
    
24.     printf("\nmatA + matB:\n");
25.     matrix_Print(matrix_Add(matA,matB));
26. 
    
27.     printf("\nmatA - matB:\n");
28.     matrix_Print(matrix_Sub(matA,matB));
29. 
    
30.     printf("\nmatA * matB:\n");
31.     matrix_Print(matrix_Mult(matA,matB));
32. 
    
33.     printf("\n5 * matA:\n");
34.     matrix_Print(matrix_ScalarMultiply(5,matA));
35. 
    
36.     printf("\nTranspose ofmatA :\n");
37.     matrix_Print(matrix_tran(matA));
38. 
    
39.     printf("\nDeterminant of matA:\t");
40.     printf("%f",matrix_Det(matA));
41. 
    
42.     printf("\nCofactor matrix of matA :\n");
43.     matrix_Print(matrix_Cofactor(matA));
44.     printf("\nInverse matrix of matA:\n");
45.     matrix_Print(matrix_Inverse(matA));
46.     printf("\nMinor matrix of ( i , j) = (1,1) of matA:\n");
47.     matrix_Print(matrix_Minor(matA,1,1));
48. 
    
49.     return 0;
50. }

        Once we wrote a library with generic functions, it easy to develop some projects easily based on our own library. Even though, there are many efficient scientific and mathematical libraries available free of charge.  We learn Programming easily through writing programs like this. 

Happy coding!!!