Lab 8
Lab 8
Lab 8
1. Polynomial Regression
Use the same least-squares procedure (Lab #7), but now using a polynomial of order : = 0 + 1 1 + 2 2 + . . . + For the cuadratic case ( = 2): = 0 2 + 1 + 2 the sum of squares of residuals is:
2 ]
= [
=1
= [ (0 2 + 1 + 2 )]
=1
where is the number of points in the data set we are trying to fit to the model. To find the coefficients {} that minimize , take its derivative with respect to each coefficient: = 2 2 [ (0 2 + 1 + 2 )] 0
=1
= 2 [ (0 2 + 1 + 2 )] 1
=1
= 2 [ (0 2 + 1 + 2 )] 2
=1
then equate them to zero, and rearrange to get a system of linear equations:
( 2 ) 0 =1
+ ( ) 1 + ()2 =
=1 =1
( 3 ) 0 + ( 2 ) 1 + ( ) 2 =
=1 =1 =1 =1
( 4 ) 0 + ( 3 ) 1 + ( 2 ) 2 = 2
=1 =1 =1 =1
* Notes based on Chapter 15 of Applied Numerical Methods with MATLAB for Engineers and Scientists, Steven C. Chapra, 3rd Edition, McGraw-Hill, 2012.
The coefficient of determination still has the same definition as in Lab #7: 2 =
= [
=1
= [ (0 + 1 1, + 2 2, )]
=1
= 2 1, [ (0 + 1 1, + 2 2, )] 1
=1
= 2 2, [ (0 + 1 1, + 2 2, )] 2
=1
( 2, ) 2 + ( 1, ) 1 + ()0 =
=1 =1 =1 2 ( 1, 2, ) 2 + ( 1, ) 1 =1 =1 2 ( 2, ) 2 + ( 1, 2, ) 1 =1 =1
+ ( 1, ) 0 = 1,
=1 =1
+ ( 2, ) 0 = 2,
=1 =1
[0,1
1,1
0,1 0,2 [ 0,
1,1 1,2 1,
In order for us to be able to solve this system of equations, the matrix [] has to be a squared matrix ( = + 1). This means that if we are fitting a parabola ( = 2) to data, we can only use = 2 + 1 = 3 data points to obtain the coefficients {}. Instead, we will minimize the sum of the squares of residuals :
2
= ( , )
=1 =0
The minimization yields: [[]T []] {} = {[]T {}} This means that we can solve for {} by multiplying on both sides of the equation by the inverse of [[]T []]: [[]T []]
1
{[]T {}}
{[]T {}}
{[]T {}}
and the estimates of y ( ) that correspond to the original data points () can be calculated using = []{} The calculation of the coefficient of determination 2 , follows the usual definition of 2 = With Python: Once we have calculated the matrix [] and the vector {}, we can solve for the coefficients {} by explicitly taking the inverse of [[]T []] and multiplying it by {[]T {}}: from numpy.linalg import inv a = inv(transpose(Z)*Z)*(transpose(Z)*y) or by using the function lstsq: from numpy.linalg import lstsq a = numpy.linalg.lstsq(Z,y,rcond=-1) Returns the least-squares solution to a linear matrix equation. Solves the equation = by computing a vector that minimizes the Euclidean 2-norm | |2 . The equation may be under-, well-, or over- determined (i.e., the number of linearly independent rows of a can be less than, equal to, or greater than its number of linearly independent columns). If is square and of full rank, then (but for round-off error) is the exact solution of the equation. Parameters: Z: (M, N) array like Coefficient matrix. y: {(M,), (M, K)} array_like Ordinate or dependent variable values. If y is two-dimensional, the leastsquares solution is calculated for each of the K columns of y. rcond: float, optional
.
Returns:
Cut-off ratio for small singular values of Z. Singular values are set to zero if they are smaller than rcond times the largest singular value of Z. a: {(N,), (N, K)} ndarray Least-squares solution. If is two-dimensional, the solutions are in the K columns of . residuals: {(), (1,), (K,)} ndarray Sums of residuals; squared Euclidean 2-norm for each column in . If the rank of is < N or > M, this is an empty array. If is 1-dimensional, this is a (1,) shape array. Otherwise the shape is (K). rank: int Rank of matrix . s: (min(M, N),) ndarray Singular values of .
4. Non-Linear Regression
First try linearizing the model (like we did in Lab #7). If model is not linearizable, perform non-linear regression to directly determine the least-squares fit. For example, lets say we are trying to fit the following equation to our data: = 0 (1 1 ) Step 1: Calculate the sum of squares of residuals
2 ]
= [
=1
= [ 0 (1 1 )]2
=1
and save it as a Python function. def sr(a,x,y,): Step 2: Use Pythonss fmin function to minimize by varying 0 and 1 . The general syntax is a = fmin(sr,a0,args=(x,y,),xtol=0.0001,ftol=0.0001, maxiter=None,maxfun=None,disp=0) where: a = vector of {} that minimizes sr = Python function containing expression for a0 = vector of initial guesses for {} xtol,ftol,maxiter,maxfun = values for optimization parameters args = x and y values for original data, plus additional parameters for function sr
5. Example
Given the following data: x y 10 25 20 70 30 380 40 550 50 610 60 1220 70 830 80 1450
Fit a power model ( = ), a) by linearizing the model and doing a general linear least squares fit, and b) by doing a non-linear fit a) general linear least squares fit Take the log10 on both sides to linearize the power model: log = log log = log + log log = log + log When comparing this equation with the general form of the linear model: = 0 0, + 1 1, + 2 2, + + , We get that: = log 0 = log 0, = 1 1 = 1, = log
With this, we now build our [] matrix and {} vector. For our 8 data points we should have: 1 1 1 1 1 1 1 [1 log 1 log 1 log 2 log 2 log 3 log 3 log 4 log 4 0 { }= log 5 1 log 5 log 6 log 6 log 7 log 7 log 8 ] {log 8 } 1 1 1 1 1 1 1 [1 log 10 log 25 log 20 log 70 log 30 log 380 log 40 0 log 550 { }= log 50 1 log 610 log 60 log 1220 log 70 log 830 log 80] {log 1450}
1
or
Now we can solve for {} by using the equation {} = [[]T []] function.
With the following script, we can implement this calculation into Python: from numpy import array,matrix,transpose,append,ones,log10,var from numpy.linalg import inv,lstsq from pylab import figure,plot,grid,xlabel,ylabel,legend # Enter data x = array([10.,20.,30.,40.,50.,60.,70.,80.]) y = array([25.,70.,380.,550.,610.,1220.,830.,1450.]) # Calculate statistics for the data n = len(y) St = var(y)*n # 1. Linearized fit using general model z0 = transpose(matrix(ones(n))) z1 = transpose(matrix(log10(x))) Z = append(z0,z1,axis=1) ylin = transpose(matrix(log10(y))) # a. Using matrix inverse a = inv(transpose(Z)*Z)*(transpose(Z)*ylin) ylin_est = Z*a y_est = 10.**(array(transpose(ylin_est))[0]) figure() plot(x,y,'o',x,y_est) alpha = 10.**a[0,0] beta = a[1,0] Sr = sum((y-y_est)**2) r2 = (St-Sr)/St print 'Fit to linearized model' print ' Using matrix inversion:' print ' alpha =',alpha print ' beta =',beta print ' r2 =',r2 # b. Using lstsq function a = matrix(lstsq(Z,ylin)[0]) ylin_est = Z*a y_est = 10.**(array(transpose(ylin_est))[0]) plot(x,y_est,'--') alpha = 10.**a[0,0] beta = a[1,0] Sr = sum((y-y_est)**2) r2 = (St-Sr)/St print ' Using lstsq function' print ' alpha =',alpha print ' beta =',beta print ' r2 =',r2
The output from this script is: Fit to linearized model Using matrix inversion: alpha = 0.274137342013 beta = 1.98417625576 r2 = 0.808818120972 Using lstsq function alpha = 0.274137342013 beta = 1.98417625576 r2 = 0.808818120972 Fitting of the linearized version of the power model gives us a coefficient of determination of 2 = 0.80882 Let us verify this result by using the polyfit tool. Add the following lines to the script above: # c. Using polyfit from numpy import polyfit,polyval xlin = log10(x) ylin = log10(y) p = polyfit(xlin,ylin,1) ylin_est = polyval(p,xlin) y_est = 10.**(array(ylin_est)) plot(x,y_est,':k') alpha = 10.**p[1] beta = p[0] Sr = sum((y-y_est)**2) r2 = (St-Sr)/St print ' Using polyfit function' print ' alpha =',alpha print ' beta =',beta print ' r2 =',r2 The output from this part is: Using polyfit function alpha = 0.274137342013 beta = 1.98417625576 r2 = 0.808818120972 This gives us the same value for the coefficient of determination of 2 = 0.80882 Notice that this linearized version of the power model gives us a worst fit than the linear fit to a straight line, that we did on this same data last week, which gave us 2 = 0.88048.
b) non-linear fit To do the non-linear fit of the power model, we first have to create a Python function containing the sum of squared residuals:
= [
=1
2 ]
= [ ]
=1
By setting 0 = and 1 = , we can create the following function: def sr(a,xm,ym): yp = a[0]*(xm**a[1]) return sum((ym-yp)**2) Now we are ready to use Pythons fmin tool to find the values of {} that minimize this quantity. Add the following lines to the script (notice that here I use the values obtained for and from the linear fit, as initial guesses for the fmin search): from scipy.optimize import fmin a0 = array([alpha,beta]) a = fmin(sr,a0,args=(x,y), \ xtol=0.0001,ftol=0.0001, \ maxiter=None,maxfun=None, \ disp=0) alpha = a[0] beta = a[1] y_est = alpha*(x**beta) plot(x,y_est) Sr = sr(a,x,y) r2 = (St-Sr)/St print 'Fit to non-linear model:' print ' Using fmin function' print ' alpha =',alpha print ' beta =',beta print ' r2 =',r2 # Add labels and legend to the plot grid() xlabel('x') ylabel('y') legend(('data','fit to linearized data using $inv$', \ 'fit to linearized data using $lstsq$', \ 'fit to linearized data using $polyfit$', \ 'fit to non-linear data'),loc='upper left') The output now shows: Fit to non-linear model: Using fmin function alpha = 2.53837245238
beta = 1.43585549642 r2 = 0.876898062742 This non-linear fit of the power model gives us a coefficient of determination of 2 = 0.8769, so it is still not as good as the plain straight-line model, but it is an improvement over the linearized version of this same model. The script finishes by plotting all the fits:
6. Problem
The following data were collected for the steady flow of water in a concrete circular pipe: diameter, (m) 0.3 0.6 0.9 0.3 0.6 0.9 0.3 0.6 0.9 slope, (m/m) 0.001 0.001 0.001 0.01 0.01 0.01 0.05 0.05 0.05 flow, (m3/s) 0.04 0.24 0.69 0.13 0.82 2.38 0.31 1.95 5.66 Fit the following model to this data: = 0 1 2 , using: a) generalized linear least squares regression and b) non-linear regression and compare the values for the parameters {} and the coefficient of determination obtained for these. Plot your fits along with the data. Steps to follow: 1. Linearize the model. 2. Perform a general linear least squares regression on the data, and calculate the coefficient of determination for this fit. 3. Perform a non-linear least squares regression on the data, and calculate the coefficient of determination for this fit. 4. Plot the data along with the fits.
The Python code used to solve this problem follows: #Fit a power model Q = a0*D^(a1)S^(a2) to the data using #linearization and non-linear regression. Compare the #coefficient of variation for these. from from from from from from # D S Q numpy import array,matrix,transpose,append numpy import ones,log10,var,linspace,meshgrid numpy.linalg import inv pylab import figure matplotlib import cm scipy.optimize import fmin
# Calculate statistics for the data n = len(Q) St = var(Q)*n #Fit to linearized data using general linear least squares Z = transpose(matrix(ones(n))) Z = append(Z,transpose(matrix(log10(D))),axis=1) Z = append(Z,transpose(matrix(log10(S))),axis=1) Qlin = transpose(matrix(log10(Q))) a = inv(transpose(Z)*Z)*(transpose(Z)*Qlin) Qlin_est = Z*a Q_est = 10.**(array(transpose(Qlin_est))[0]) a0l = 10.**a[0,0] a1l = a[1,0] a2l = a[2,0] Sr = sum((Q-Q_est)**2) r2 = (St-Sr)/St print 'Fit to linearized model' print ' a0 =',a0l print ' a1 =',a1l print ' a2 =',a2l print ' r2 =',r2 # Non-linear fit def sr(a,dm,sm,qm): qp = a[0]*(dm**a[1])*(sm**a[2]) return sum((qm-qp)**2) a0 = array([a0l,a1l,a2l]) a = fmin(sr,a0,args=(D,S,Q), \ xtol=0.0001,ftol=0.0001, \ maxiter=None,maxfun=None, \ disp=0) a0nl = a[0]
a1nl = a[1] a2nl = a[2] Sr = sr(a,D,S,Q) r2 = (St-Sr)/St print 'Fit to non-linear model:' print ' a0 =',a0nl print ' a1 =',a1nl print ' a2 =',a2nl print ' r2 =',r2 # Plot fits with data x = linspace(0.2,1.0,40) y = linspace(0.0,0.051,40) [X,Y] = meshgrid(x,y) Z1 = a0l*(X**a1l)*(Y**a2l) Z2 = a0nl*(X**a1nl)*(Y**a2nl) fig = figure(figsize=(14,7)) ax1 = fig.add_subplot(1,2,1,projection='3d') ax1.plot_surface(X,Y,Z1,rstride=1,cstride=1, \ cmap=cm.coolwarm,linewidth=0, \ antialiased=False) ax1.set_zlim(0.0,6.0) ax1.scatter(D,S,Q) ax1.set_xlabel('D (m)') ax1.set_ylabel('S (m/m)') ax1.set_zlabel('Q (m$^{3}$/s)') ax2 = fig.add_subplot(1,2,2,projection='3d') ax2.plot_surface(X,Y,Z2,rstride=1,cstride=1, \ cmap=cm.coolwarm,linewidth=0, \ antialiased=False) ax2.set_zlim(0.0,6.0) ax2.scatter(D,S,Q) ax2.set_xlabel('D (m)') ax2.set_ylabel('S (m/m)') ax2.set_zlabel('Q (m$^{3}$/s)') The parameters and 2 values resulting from the fits are: Fit to a0 = a1 = a2 = r2 = Fit to a0 = a1 = a2 = r2 = linearized model 36.3813323425 2.62793693821 0.531987421319 0.999862616234 non-linear model: 37.4297395804 2.63046170577 0.538054665112 0.999998444143