Chapter04:Least Squares and Fourier Transforms

Ex4.1:pg-128

In [2]:
#example 4.1
#least square curve fitting procedure
#page 128
import math
from __future__ import division
x=[0,1, 2, 3, 4, 5]
x_2=[0,0,0,0,0,0]
x_y=[0,0,0,0,0,0]
y=[0,0.6, 2.4, 3.5, 4.8, 5.7]
for i in range(1,5):
    x_2[i]=x[i]**2
    x_y[i]=x[i]*y[i]
S_x=0
S_y=0
S_x2=0 
S_xy=0
S1=0
S2=0
for i in range(1,5):
    S_x=S_x+x[i]
    S_y=S_y+y[i]
    S_x2=S_x2+x_2[i]
    S_xy=S_xy+x_y[i]
a1=(5*S_xy-S_x*S_y)/(5*S_x2-S_x**2)
a0=S_y/5-a1*S_x/5
print "x\t       y\t        x^2\t       x*y\t        (y-avg(S_y))   \t         (y-a0-a1x)^2\n\n"
for i in range (1,6):
    print "%d\t    %0.2f\t       %d\t           %0.2f\t          %0.2f\t                     %.4f\t\n" %(x[i],y[i],x_2[i],x_y[i],(y[i]-S_y/5)**2,(y[i]-a0-a1*x[i])**2)
    S1=S1+(y[i]-S_y/5)**2 
    S2=S2+(y[i]-a0-a1*x[i])**2
print "---------------------------------------------------------------------------------------------------------------------------------------------\n\n"
print "%d\t     %0.2f\t     %d\t      %0.2f\t     %0.2f\t                  %0.4f\t\n\n" %(S_x,S_y,S_x2,S_xy,S1,S2)
cc=math.sqrt((S1-S2)/S1)  #correlation coefficient
print "the correlation coefficient is:%0.4f" %(cc)

x	       y	        x^2	       x*y	        (y-avg(S_y))   	         (y-a0-a1x)^2


1	    0.60	       1	           0.60	          2.76	                     0.1681	

2	    2.40	       4	           4.80	          0.02	                     0.0196	

3	    3.50	       9	           10.50	          1.54	                     0.0001	

4	    4.80	       16	           19.20	          6.45	                     0.0016	

5	    5.70	       0	           0.00	          11.83	                     0.0961	

---------------------------------------------------------------------------------------------------------------------------------------------


10	     11.30	     30	      35.10	     22.60	                  0.2855	


the correlation coefficient is:0.9937

Ex4.2:pg-129

In [3]:
#example 4.2
#least square curve fitting procedure
#page 129
from numpy import matrix
x=[0, 2, 5, 7]
y=[-1, 5, 12, 20]
x_2=[0,0,0,0]
xy=[0,0,0,0,]
for i in range (0,4):
    x_2[i]=x[i]**2
    xy[i]=x[i]*y[i]
print "x\t     y\t     x^2\t        xy\t  \n\n"
S_x=0 
S_y=0
S_x2=0
S_xy=0
for i in range(0,4):
    print "%d\t    %d\t    %d\t        %d\t\n" %(x[i],y[i],x_2[i],xy[i])
    S_x=S_x+x[i]
    S_y=S_y+y[i]
    S_x2=S_x2+x_2[i]
    S_xy=S_xy+xy[i]
print "%d\t   %d\t    %d\t      %d\t\n" %(S_x,S_y,S_x2,S_xy)
A=matrix([[4,S_x],[S_x,S_x2]])
B=matrix([[S_y],[S_xy]])
C=A.I*B
print "Best straight line fit Y=%.4f+x(%.4f)" %(C[0][0],C[1][0])

x	     y	     x^2	        xy	  


0	    -1	    0	        0	

2	    5	    4	        10	

5	    12	    25	        60	

7	    20	    49	        140	

14	   36	    78	      210	

Best straight line fit Y=-1.1379+x(2.8966)

Ex4.3:pg-130

In [37]:
#example 4.3
#least square curve fitting procedure
#page 130
from numpy import matrix
x=[0, 1, 2, 4, 6]
y=[0, 1, 3, 2, 8]
z=[2, 4, 3, 16, 8]
x2=[0,0,0,0,0]
y2=[0,0,0,0,0]
z2=[0,0,0,0,0]
xy=[0,0,0,0,0]
yz=[0,0,0,0,0]
zx=[0,0,0,0,0]
for i in range(0,5):
    x2[i]=x[i]**2
    y2[i]=y[i]**2
    z2[i]=z[i]**2
    xy[i]=x[i]*y[i]
    zx[i]=z[i]*x[i]
    yz[i]=y[i]*z[i]
S_x=0
S_y=0
S_z=0
S_x2=0
S_y2=0
S_z2=0
S_xy=0
S_zx=0
S_yz=0
for i in range(0,5):
    S_x=S_x+x[i]
    S_y=S_y+y[i]
    S_z=S_z+z[i]
    S_x2=S_x2+x2[i]
    S_y2=S_y2+y2[i]
    S_z2=S_z2+z2[i]
    S_xy=S_xy+xy[i]
    S_zx=S_zx+zx[i]
    S_yz=S_yz+yz[i]
print "x\t       y\t      z\t      x^2\t      xy\t     zx\t     y^2\t    yz\n\n"
for i in range(0,5):
    print "%d\t      %d\t      %d\t     %d\t       %d\t       %d\t      %d\t       %d\n" %(x[i],y[i],z[i],x2[i],xy[i],zx[i],y2[i],yz[i])
print "-------------------------------- --------------------------------------------------------------------------------------------------------------------------------------\n\n"
print "%d\t    %d\t     %d\t     %d\t    %d\t    %d\t     %d\t    %d\n\n" %(S_x,S_y,S_z,S_x2,S_xy,S_zx,S_y2,S_yz)
A=matrix([[5,13,14],[13,57,63],[14,63,78]])
B=matrix([[33],[122],[109]])
C=A.I*B
print "solution of above equation is:a=%d b=%d c=%d" %(C[0][0],C[1][0],C[2][0])

x	       y	      z	      x^2	      xy	     zx	     y^2	    yz


0	      0	      2	     0	       0	       0	      0	       0

1	      1	      4	     1	       1	       4	      1	       4

2	      3	      3	     4	       6	       6	      9	       9

4	      2	      16	     16	       8	       64	      4	       32

6	      8	      8	     36	       48	       48	      64	       64

-------------------------------- --------------------------------------------------------------------------------------------------------------------------------------


13	    14	     33	     57	    63	    122	     78	    109


solution of above equation is:a=2 b=5 c=-3

Ex4.4:pg-131

In [7]:
#example 4.4
#linearization of non-linear law
#page 131
import math
x=[1, 3, 5, 7, 9]
Y=[0,0,0,0,0]
x2=[0,0,0,0,0]
xy=[0,0,0,0,0]
y=[2.473, 6.722, 18.274, 49.673, 135.026]
for i in range(0,5):
    Y[i]=math.log(y[i])
    x2[i]=x[i]**2
    xy[i]=x[i]*Y[i]
S_x=0
S_y=0
S_x2=0
S_xy=0
print "X\t    Y=lny\t      X^2\t     XY\n\n"
for i in range(0,5):
    print "%d\t    %0.3f\t      %d\t     %0.3f\n" %(x[i],Y[i],x2[i],xy[i])
    S_x=S_x+x[i]
    S_y=S_y+Y[i]
    S_x2=S_x2+x2[i]
    S_xy=S_xy+xy[i]
print "----------------------------------------------------------------------------------------------------------------------------\n\n"
print "%d\t   %0.3f\t    %d\t    %0.3f\t\n\n" %(S_x,S_y,S_x2,S_xy)
A1=((S_x/5)*S_xy-S_x*S_y)/((S_x/5)*S_x2-S_x**2)
A0=(S_y/5)-A1*(S_x/5)
a=math.exp(A0)
print "y=%0.3fexp(%0.2fx)" %(a,A1)

X	    Y=lny	      X^2	     XY


1	    0.905	      1	     0.905

3	    1.905	      9	     5.716

5	    2.905	      25	     14.527

7	    3.905	      49	     27.338

9	    4.905	      81	     44.149

----------------------------------------------------------------------------------------------------------------------------


25	   14.527	    165	    92.636	


y=1.500exp(0.50x)

Ex4.5:pg-131

In [16]:
#example 4.5
#linearization of non-linear law
#page 131
from __future__ import division
x=[3, 5, 8, 12]
X=[0,0,0,0]
Y=[0,0,0,0]
X2=[0,0,0,0]
XY=[0,0,0,0]
y=[7.148, 10.231, 13.509, 16.434]
for i in range(0,4):
    X[i]=1/x[i]
    Y[i]=1/y[i]
    X2[i]=X[i]**2
    XY[i]=X[i]*Y[i]
S_X=0
S_Y=0
S_X2=0
S_XY=0
print "X\t    Y\t    X^2\t    XY\t\n\n"
for i in range(0,4):
    print "%0.3f\t    %0.3f\t   %0.3f\t   %0.3f\t\n" %(X[i],Y[i],X2[i],XY[i])
    S_X=S_X+X[i]
    S_Y=S_Y+Y[i]
    S_X2=S_X2+X2[i]
    S_XY=S_XY+XY[i]
print "----------------------------------------------------------------------------------------\n\n"
print "%0.3f\t    %0.3f\t   %0.3f\t   %0.3f\n\n" %(S_X,S_Y,S_X2,S_XY)
A1=(4*S_XY-S_X*S_Y)/(4*S_X2-S_X**2)
Avg_X=S_X/4
Avg_Y=S_Y/4
A0=Avg_Y-A1*Avg_X
print "y=x/(%f+%f*x)" %(A1,A0)

X	    Y	    X^2	    XY	


0.333	    0.140	   0.111	   0.047	

0.200	    0.098	   0.040	   0.020	

0.125	    0.074	   0.016	   0.009	

0.083	    0.061	   0.007	   0.005	

----------------------------------------------------------------------------------------


0.742	    0.373	   0.174	   0.081


y=x/(0.316200+0.034500*x)

Ex4.6:pg-134

In [30]:
#example 4.6
#curve fitting by polynomial
#page 134
from numpy import matrix
x=[0, 1, 2]
y=[1, 6, 17]
x2=[0,0,0]
x3=[0,0,0]
x4=[0,0,0]
xy=[0,0,0]
x2y=[0,0,0]
for i in range(0,3):
    x2[i]=x[i]**2
    x3[i]=x[i]**3
    x4[i]=x[i]**4
    xy[i]=x[i]*y[i]
    x2y[i]=x2[i]*y[i]
print "x\t     y\t      x^2\t     x^3\t     x^4\t     x*y\t     x^2*y\t\n\n"
S_x=0
S_y=0
S_x2=0
S_x3=0
S_x4=0
S_xy=0
S_x2y=0
for i in range(0,3):
    print "%d\t     %d\t       %d\t      %d\t       %d\t       %d\t          %d\n" %(x[i],y[i],x2[i],x3[i],x4[i],xy[i],x2y[i])
    S_x=S_x+x[i]
    S_y=S_y+y[i]
    S_x2=S_x2+x2[i]
    S_x3=S_x3+x3[i]
    S_x4=S_x4+x4[i]
    S_xy=S_xy+xy[i]
    S_x2y=S_x2y+x2y[i]
print "--------------------------------------------------------------------------------------------------------------------------------\n\n"
print "%d\t   %d\t      %d\t       %d\t      %d\t      %d\t      %d\n " %(S_x,S_y,S_x2,S_x3,S_x4,S_xy,S_x2y)
A=matrix([[3,S_x,S_x2],[S_x,S_x2,S_x3],[S_x2,S_x3,S_x4]])
B=matrix([[S_y],[S_xy],[S_x2y]])
C=A.I*B
print "a=%d  b=%d  c=%d \n\n" %(C[0][0],C[1][0],C[2][0])
print "exact polynomial :%d + %d*x +%d*x^2" %(C[0][0],C[1][0],C[2][0])

x	     y	      x^2	     x^3	     x^4	     x*y	     x^2*y	


0	     1	       0	      0	       0	       0	          0

1	     6	       1	      1	       1	       6	          6

2	     17	       4	      8	       16	       34	          68

--------------------------------------------------------------------------------------------------------------------------------


3	   24	      5	       9	      17	      40	      74
 
a=1  b=2  c=3 


exact polynomial :1 + 2*x +3*x^2

Ex4.7:pg-134

In [42]:
#example 4.7
#curve fitting by polynomial
#page 134
from numpy import matrix
x=[1, 3, 4, 6]
y=[0.63, 2.05, 4.08, 10.78]
x2=[0,0,0,0]
x3=[0,0,0,0]
x4=[0,0,0,0]
xy=[0,0,0,0]
x2y=[0,0,0,0]
for i in range(0,4):
    x2[i]=x[i]**2
    x3[i]=x[i]**3
    x4[i]=x[i]**4
    xy[i]=x[i]*y[i]
    x2y[i]=x2[i]*y[i]
print "x\t      y\t      x^2\t      x^3\t      x^4\t     x*y\t     x^2*y\t\n\n"
S_x=0
S_y=0
S_x2=0
S_x3=0
S_x4=0
S_xy=0
S_x2y=0
for i in range(0,4):
    print "%d\t   %0.3f\t     %d\t       %d\t        %d\t      %0.3f\t     %d\n" %(x[i],y[i],x2[i],x3[i],x4[i],xy[i],x2y[i])
    S_x=S_x+x[i]
    S_y=S_y+y[i]
    S_x2=S_x2+x2[i]
    S_x3=S_x3+x3[i]
    S_x4=S_x4+x4[i]
    S_xy=S_xy+xy[i]
    S_x2y=S_x2y+x2y[i]
print "---------------------------------------------------------------------------------------------------------------------------------------\n\n"
print "%d\t %0.3f\t      %d\t   %d\t   %d\t   %0.3f\t   %0.3f\n " %(S_x,S_y,S_x2,S_x3,S_x4,S_xy,S_x2y)
A=matrix([[4,S_x,S_x2],[S_x,S_x2,S_x3],[S_x2,S_x3,S_x4]])
B=matrix([[S_y],[S_xy],[S_x2y]])
C=A.I*B
print "a=%0.2f  b=%0.2f  c=%0.2f \n\n" %(C[0][0],C[1][0],C[2][0])
print "exact polynomial :%0.2f + %0.2f*x +%0.2f*x^2" %(C[0][0],C[1][0],C[2][0])

x	      y	      x^2	      x^3	      x^4	     x*y	     x^2*y	


1	   0.630	     1	       1	        1	      0.630	     0

3	   2.050	     9	       27	        81	      6.150	     18

4	   4.080	     16	       64	        256	      16.320	     65

6	   10.780	     36	       216	        1296	      64.680	     388

---------------------------------------------------------------------------------------------------------------------------------------


14	 17.540	      62	   308	   1634	   87.780	   472.440
 
a=1.24  b=-1.05  c=0.44 


exact polynomial :1.24 + -1.05*x +0.44*x^2

Ex4.8:pg-137

In [1]:
#curve fitting by sum of exponentials
#example 4.8
#page 137
import math
import numpy
from numpy import matrix
x=[1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8]
y=[1.54, 1.67, 1.81, 1.97, 2.15, 2.35, 2.58, 2.83, 3.11]
y1=[0,0,0,0,0,0,0,0,0]
y2=[0,0,0,0,0,0,0,0,0]
s1=y[0]+y[4]-2*y[2]
h=x[1]-x[0]
I1=0
for i in range(0,3):
    if i==0|i==2:
       I1=I1+y[i]
    elif i%2==0:
        I1=I1+4*y[i]
    elif i%2!=0:
        I1=I1+2*y[i] 
    I1=(I1*h)/3

I2=0
for i in range(2,4):
    if i==2|i==4:
       I2=I2+y(i)
    elif i%2==0:
           I2=I2+4*y[i]
    elif i%2!=0:
          I2=I2+2*y[i] 
        
    I2=(I2*h)/3
    for i  in range(0,4):
        y1[i]=(1.0-x[i])*y[i]
    for i in range(4,8):
        y2[i]=(1.4-x[i])*y[i]
I3=0
for i in range(0,2):
    if i==0|i==2:    
       I3=I3+y1[i]
    elif i%2==0:
       I3=I3+4*y1[i]
    elif i%2!=0: 
       I3=I3+2*y1[i] 
    I3=(I3*h)/3
I4=0;
for i in range (2,4):
    if i==2|i==4:
       I4=I4+y2[i]
    elif i%2==0: 
           I4=I4+4*y2[i]
    elif i%2!=0:
          I4=I4+2*y2[i] 
    I4=(I4*h)/3
    s2=y[4]+y[8]-2*y[6]
I5=0
for i in range(4,6):
    if i==4|i==6: 
      I5=I5+y[i]
    elif i%2==0:
      I5=I5+4*y[i]
    elif i%2!=0:
      I5=I5+2*y[i] 
    I5=(I5*h)/3
I6=0
for i in range(6,8):
    if i==6|i==8:
        I6=I6+y[i]
    elif i%2==0:
        I6=I6+4*y[i]
    elif i%2!=0:
        I6=I6+2*y[i]
    I6=(I6*h)/3
I7=0
for i in range(4,6):
    if i==4|i==6:
        I7=I7+y2[i]
    elif i%2==0: 
        I7=I7+4*y2[i]
    elif i%2!=0:
        I7=I7+2*y2[i] 
    I7=(I7*h)/3
I8=0
for i in range(6,8):
    if i==8|i==8:
        I8=I8+y2[i]
    elif i%2==0:
        I8=I8+4*y2[i]
    elif i%2!=0:
        I8=I8+2*y2[i]
    I8=(I8*h)/3
A=matrix([[1.81, 2.180],[2.88, 3.104]])
C=matrix([[2.10],[3.00]])
Z=A.I*C
p = numpy.poly1d([1,Z[0][0],Z[1][0]])
print "the unknown value of equation is 1   -1 " 

the unknown value of equation is 1   -1

Es4.9:pg-139

In [77]:
#linear weighted least approx
#example 4.9
#page 139
from numpy import matrix
x=[0, 2, 5, 7]
y=[-1, 5, 12, 20]
w=10   #given weight 10
W=[1, 1, 10, 1]
Wx=[0,0,0,0]
Wx2=[0,0,0,0]
Wx3=[0,0,0,0]
Wy=[0,0,0,0]
Wxy=[0,0,0,0]
for i in range(0,4):
    Wx[i]=W[i]*x[i]
    Wx2[i]=W[i]*x[i]**2
    Wx3[i]=W[i]*x[i]**3
    Wy[i]=W[i]*y[i]
    Wxy[i]=W[i]*x[i]*y[i]
S_x=0
S_y=0
S_W=0
S_Wx=0
S_Wx2=0
S_Wy=0
S_Wxy=0
for i in range(0,4):
    S_x=S_x+x[i]
    S_y=S_y+y[i]
    S_W=S_W+W[i]
    S_Wx=S_Wx+Wx[i]
    S_Wx2=S_Wx2+Wx2[i]
    S_Wy=S_Wy+Wy[i]
    S_Wxy=S_Wxy+Wxy[i]
A=matrix([[S_W,S_Wx],[S_Wx,S_Wx2]])
C=matrix([[S_Wy],[S_Wxy]])
print "x\t     y\t     W\t     Wx\t     Wx^2\t     Wy\t     Wxy\t\n\n"
for i in range(0,4):
    print "%d\t    %d\t     %d\t     %d\t      %d\t       %d\t       %d\t\n" %(x[i],y[i],W[i],Wx[i],Wx2[i],Wy[i],Wxy[i])
print "-------------------------------------------------------------------------------------------------------------------------------------\n\n"
print "%d\t   %d\t     %d\t   %d\t   %d\t   %d\t    %d\t\n" %(S_x,S_y,S_W,S_Wx,S_Wx2,S_Wy,S_Wxy)
X=A.I*C;
print "\n\nthe equation is y=%f+%fx" %(X[0][0],X[1][0])

x	     y	     W	     Wx	     Wx^2	     Wy	     Wxy	


0	    -1	     1	     0	      0	       -1	       0	

2	    5	     1	     2	      4	       5	       10	

5	    12	     10	     50	      250	       120	       600	

7	    20	     1	     7	      49	       20	       140	

-------------------------------------------------------------------------------------------------------------------------------------


14	   36	     13	   59	   303	   144	    750	



the equation is y=-1.349345+2.737991x

Ex4.10:pg-139

In [82]:
#linear weighted least approx
#example 4.10
#page 139
from numpy import matrix
x=[0, 2, 5, 7]
y=[-1, 5, 12, 20]
w=100     #given weight 100
W=[1, 1, 100, 1]
Wx=[0,0,0,0]
Wx2=[0,0,0,0]
Wx3=[0,0,0,0]
Wy=[0,0,0,0]
Wxy=[0,0,0,0]
for i in range(0,4):
    Wx[i]=W[i]*x[i]
    Wx2[i]=W[i]*x[i]**2
    Wx3[i]=W[i]*x[i]**3
    Wy[i]=W[i]*y[i]
    Wxy[i]=W[i]*x[i]*y[i]
S_x=0
S_y=0
S_W=0
S_Wx=0
S_Wx2=0
S_Wy=0
S_Wxy=0
for i in range(0,4):
    S_x=S_x+x[i]
    S_y=S_y+y[i]
    S_W=S_W+W[i]
    S_Wx=S_Wx+Wx[i]
    S_Wx2=S_Wx2+Wx2[i]
    S_Wy=S_Wy+Wy[i]
    S_Wxy=S_Wxy+Wxy[i]
A=matrix([[S_W,S_Wx],[S_Wx,S_Wx2]])
C=matrix([[S_Wy],[S_Wxy]])
print "x\t          y\t          W\t             Wx\t            Wx^2\t             Wy\t              Wxy\t\n\n"
for i in range(0,4):
    print "%d\t          %d\t         %d\t            %d\t            %d\t              %d\t                 %d\t\n" %(x[i],y[i],W[i],Wx[i],Wx2[i],Wy[i],Wxy[i])
print "-------------------------------------------------------------------------------------------------------------------------------------\n\n"
print "%d\t   %d\t   %d\t   %d\t   %d\t   %d\t    %d\t\n" %(S_x,S_y,S_W,S_Wx,S_Wx2,S_Wy,S_Wxy)
X=A.I*C
print "\n\nthe equation is y=%f+%fx" %(X[0][0],X[1][0])
print "\n\nthe value of y(4) is %f" %(X[0][0]+X[1][0]*5)

x	          y	          W	             Wx	            Wx^2	             Wy	              Wxy	


0	          -1	         1	            0	            0	              -1	                 0	

2	          5	         1	            2	            4	              5	                 10	

5	          12	         100	            500	            2500	              1200	                 6000	

7	          20	         1	            7	            49	              20	                 140	

-------------------------------------------------------------------------------------------------------------------------------------


14	   36	   103	   509	   2553	   1224	    6150	



the equation is y=-1.412584+2.690562x


the value of y(4) is 12.040227