Power Point to accompany Introduction to MATLAB 7

Agenda Linear Least Squares Curve Fitting Using polyfit for polynomial LLS Function Discovery Linear, Power, or Exponential? Adapting polyfit to nonstandard curve fitting Regression Analysis 3 D Plotting

Data Modeling from http: //www. cs. princeton. edu/courses/archive/fall 09/cos 323/lectures/cos 323_s 06_lecture 0 7_lsq. ppt • Given: data points, functional form, find constants in function • Example: given (xi, yi), find line through them; i. e. , find a and b in y = ax+b (x 6, y 6) (x 3, y 3) y=ax+b (x 5, y 5) (x 1, y 1) (x 7, y 7) (x 2, y 2) (x 4, y 4)

Data Modeling • You might do this because you actually care about those numbers… – Example: measure position of falling object, fit parabola position time p = – 1/2 gt 2 Estimate g from fit

Data Modeling • … or because some aspect of behavior is unknown and you want to ignore it – Example: measuring relative resonant frequency of two ions, want to ignore magnetic field drift

Least Squares • Nearly universal formulation of fitting: minimize squares of differences between data and function – Example: for fitting a line, minimize with respect to a and b – Most general solution technique: take derivatives w. r. t. unknown variables, set equal to zero – Another Approach: use linear algebra • solve overconstrained system of equations

One Scenario We measured the following data in a physics lab on Hooke's Law (spring force) Stretch Distance, m Force, N 0 0 0. 1 10. 1 0. 25 24. 5 0. 3 32 0. 6 55 We know F = kx and want to estimate k from our data. (Note the "noise" in data)

How to find our measured k? We can fit a straight line to the data y = mx + b or rather kx + b = F m is the slope of the line b is the y-intercept (vertical offset) Plug each of our data points into formula: x, m Force, N k* 0 +b= 0 k* 0. 1 +b= 10. 1 k* 0. 25 +b= 24. 5 k* 0. 3 +b= 32 k* 0. 6 +b= 55

Linear Algebra Form 0 1 0. 1 1 0. 25 1 0. 3 1 32 0. 6 1 55 * k b = 0 10. 1 24. 5 OR, A * u = F A is matrix of x (spring length) and "ones" u is vector of unknowns u(1)=k, u(2)=b F is vector of Force measurements

An "Overconstrained" System More equations than we have unknowns no "exact" solution but we can solve a modified system: Au=y A TA u = A Ty (ATA)-1 (ATA) u = (ATA)-1 ATy In Matlab: u = inv( A' * A) * A' * y (no '. ' we want matrix mult)

Matlab Solution therefore: k = 92. 0714 b = 1. 3021

Shortcut using polyfit >> u = polyfit(x, F, 1) u= 92. 0714 1. 3021 Again, k = u(1) = 92. 0714 b = u(2) = 1. 3021

Lets plot our data and model

polyfit, polyval p = polyfit( x, y, n); returns n'th order polynomial that's the closest fit to x-y data (x and y are vectors, n a scalar) yp = polyval( p, x ) returns the vector yp = the values of polynomial p at each value in x

Things to try Based on your model, what Force would you expect to see if the spring were stretched 0. 5 m ? What if we fit a quadratic to our data instead of a polynomial? Would that change k?

Example Curves polyfit can model Linear Quadratic Cuartic Cubic, etc

Other types of Models Exponential Power Law How can we determine these models?

Function Discovery Sometimes we don't know the underlying laws We've got a batch of experimental data contains noise perhaps irregular sampling To be useful, data must be fit by math model reduces noise provides estimates between samples How to choose math model to fit? linear? polynomial? power? exponential?

3 Common Fitting Fns The linear function y mx b The power function y bxm Its slope is m and its intercept is b. x raised to some power m m does not need to be an integer or positive no other terms – NOT a polynomial The exponential function y b(10)mx can also be written as y bemx base (10 or e) raised to some power mx In all 3, fitting involves estimating m and b !

Function Discovery – Which is it? ? ? involves plotting data on linear, semilogy and loglog scales and looking for straight line behavior The power function y 2 x 0 5 and the exponential function y 101 x Figure 5. 3– 8 5 -55

Plots tell which model to use Each function gives a straight line when plotted using a specific set of axes: 1. The linear function y mx b gives a straight line when plotted on rectilinear axes. Its slope is m and its intercept is b. 2. The power function y bxm gives a straight line when plotted on log-log axes. 3. The exponential function y b(10)mx and its equivalent form y bemx give a straight line when plotted on a semilog plot whose y-axis is logarithmic. More? See pages 299 -300. 5 -56

Steps for Function Discovery 1) Examine the data near the origin. The exponential function can never pass through the origin (unless of course b 0, which is a trivial case). (See Figure 5. 5– 1 for examples with b 1. ) The linear function can pass through the origin only if b 0. The power function can pass through the origin but only if m 0. (See Figure 5. 5– 2 for examples with b 1. ) 5 -57

Examples of exponential functions. Figure 5. 5– 1 5 -58

Examples of power functions. Figure 5. 5– 2 5 -59

Steps for Function Discovery (continued) 2. Plot the data using rectilinear scales. If it forms a straight line, then it can be represented by the linear function and you are finished. Otherwise, if you have data at x 0, then a. If y(0) 0, try the power function. b. If y(0) 0, try the exponential function. If data is not given for x 0, proceed to step 3. (continued…) 5 -60

Steps for Function Discovery (continued) 3. If you suspect a power function, plot the data using loglog scales. Only a power function will form a straight line on a log-log plot. If you suspect an exponential function, plot the data using the semilog scales. Only an exponential function will form a straight line on a semilog plot. (continued…) 5 -61

Steps for Function Discovery (continued) 4. In function discovery applications, we use the loglog and semilog plots only to identify the function type, but not to find the coefficients b and m. The reason is that it is difficult to interpolate on log scales. 5 -62

Plot and determine which curve fits best Read on to learn how to fit. . .

Fitting the Data – novel use of polyfit • The polyfit command fits polynomials • more complicated than the fits we want • We can use polyfit to solve for • linear, exponential, and power fits Command Description p = polyfit(x, y, n) Fits a polynomial of degree n to data described by the vectors x and y, where x is the independent variable. Returns a row vector p of length n + 1 that contains the polynomial coefficients in order of descending powers. 5 -63

Using Polyfit for Linear Curve Fit Syntax: p = polyfit(x, y, 1) where x and y contain the data, 1 means the polynomial to be fitted is linear, and p is the vector of coefficients, [p 1 p 2]. The linear fit results: given the model y = mx + b, the fit results indicate a best linear fit is made with y' = p 1 * x + p 2, where y' is the model's version of y. In other words, the first element p 1 of the vector p will be m, and the second element p 2 will be b. 5 -64

Fit for Problem 5 -30 a) best described by a linear function x=[25 30 35 40 45]; y=[5 260 480 745 1100]; pfit=polyfit(x, y, 1); xmodel=25: 45; % derive fit % construct fit curve ymodel=polyval(p, xmod); % y values for model plot(x, y, '*'); hold on; plot(xmodel, ymodel), xlabel('x'), ylabel('y') text(35, 200, 'p=53. 5 -1354. 5')

Plot for Ex 5 -30 a)

Using Polyfit for Power Curve Fit The power function: y bxm. In this case log 10 y m log 10 x log 10 b which has the form (5. 5– 5) w p 1 z p 2 where the polynomial variables w and z are related to the original data variables x and y by w log 10 y and z log 10 x. Thus we can find the power function that fits the data by typing p = polyfit(log 10(x), log 10(y), 1) The first element p 1 of the vector p will be m, and the second element p 2 will be log 10 b. We can find b from b 10 p 2. 5 -65

Using Polyfit for Exponential Curve Fit The exponential function: y b(10)mx. In this case log 10 y mx log 10 b which has the form (5. 5– 6) w p 1 z p 2 where the polynomial variables w and z are related to the original data variables x and y by w log 10 y and z x. We can find the exponential function that fits the data by typing p = polyfit(x, log 10(y), 1) The first element p 1 of the vector p will be m, and the second element p 2 will be log 10 b. We can find b from b 10 p 2. 5 -66 More? See pages 302 -303.

Least Squares Fitting The Least Squares Criterion: used to fit a function f (x). It minimizes the sum of the squares of the residuals, J. J m is defined as (5. 6– 1) J [ f (xi ) – yi ]2 i 1 We can use this criterion to compare the quality of the curve fit for two or more functions used to describe the same data. The function that gives the smallest J value gives the best fit. 5 -72

Illustration of the least squares criterion. Figure 5. 6– 1 5 -73

The least squares fit for the example data. Figure 5. 6– 2 5 -74 See pages 312 -315.

The polyfit function is based on the least-squares method. Its syntax is p = polyfit(x, y, n) Fits a polynomial of degree n to data described by the vectors x and y, where x is the independent variable. Returns a row vector p of length n+1 that contains the polynomial coefficients in order of descending powers. See page 315, Table 5. 6 -1. 5 -75

Regression using polynomials of first through fourth degree. Figure 5. 6– 3 The program is on pages 315 to 316. 5 -76

Beware of using polynomials of high degree. An example of a fifth-degree polynomial that passes through all six data points but exhibits large excursions between points. Figure 5. 6– 4 5 -77

Assessing the Quality of a Curve Fit: Denote the sum of the squares of the deviation of the y values from their mean y by S, which can be computed from m S i 1 5 -78 (yt – y )2 (5. 6– 2)

This formula can be used to compute another measure of the quality of the curve fit, the coefficient of determination, also known as the r-squared value. It is defined as r 2 1 J S (5. 6– 3) The value of S indicates how much the data is spread around the mean, and the value of J indicates how much of the data spread is unaccounted for by the model. Thus the ratio J / S indicates the fractional variation unaccounted for by the model. 5 -79

For a perfect fit, J 0 and thus r 2 1. Thus the closer r 2 is to 1, the better the fit. The largest r 2 can be is 1. It is possible for J to be larger than S, and thus it is possible for r 2 to be negative. Such cases, however, are indicative of a very poor model that should not be used. As a rule of thumb, a good fit accounts for at least 99 percent of the data variation. This value corresponds to r 2 0. 99. More? See pages 319 -320. 5 -80

Calculating J, S, r Suppose we have a vector of x and y data pfit = polyfit(x, y, 2) ; % a quadratic fit ymod = polyval(pfit, x) ; % the model y values error = ymod - y; J = sum(error. ^2) % sum of square of err ymean = mean(y); S = sum( (y-ymean). ^2) r 2 = 1 - J / S

Three-Dimensional Line Plots: The following program uses the plot 3 function to generate the spiral curve shown in Figure 5. 8– 1. t = [0: pi/50: 10*pi]; plot 3(exp(-0. 05*t). *sin(t), exp(0. 05*t). *cos(t), t) xlabel('x'), ylabel('y'), zlabel('z'), grid See the next slide. 5 -89

The curve x e 0. 05 t sin t, y e 0. 05 t cos t, z t plotted with the plot 3 function. Figure 5. 8– 1 More? See pages 334335. 5 -90

Surface Plots: The following session shows how to generate the surface plot of the function z xe [(x y 2)2 y 2], for 2 x 2 and 2 y 2, with a spacing of 0. 1. This plot appears in Figure 5. 8– 2. [X, Y] = meshgrid(-2: 0. 1: 2); Z = X. *exp(-((X-Y. ^2). ^2+Y. ^2)); mesh(X, Y, Z), xlabel('x'), ylabel('y'), zlabel(’z’) See the next slide. 5 -91

A plot of the surface z xe [(x y 2)2 y 2] created with the mesh function. Figure 5. 8– 2 5 -92 More? See pages 335 -336.

Contour Plots: The following session generates the contour plot of the function whose surface plot is shown in Figure 5. 8– 2; namely, z xe [(x y 2)2 y 2], for 2 x 2 and 2 y 2, with a spacing of 0. 1. This plot appears in Figure 5. 8– 3. >>[X, Y] = meshgrid(-2: 0. 1: 2); >>Z = X. *exp(-((X- Y. ^2). ^2+Y. ^2)); >>contour(X, Y, Z), xlabel('x'), ylabel('y') See the next slide. 5 -93

A contour plot of the surface z xe [(x y contour function. Figure 5. 8– 3 2)2 y 2] created with the More? See page 337. 5 -94

Three-dimensional plotting functions. Table 5. 8– 1 Function contour(x, y, z) Description Creates a contour plot. mesh(x, y, z) Creates a 3 D mesh surface plot. meshc(x, y, z) Same as mesh but draws contours under the surface. meshz(x, y, z) Same as mesh but draws vertical reference lines under the surface. surf(x, y, z) Creates a shaded 3 D mesh surface plot. surfc(x, y, z) Same as surf but draws contours under the surface. [X, Y] = meshgrid(x, y) Creates the matrices X and Y from the vectors x and y to define a rectangular grid. [X, Y] = meshgrid(x) Same as [X, Y]= meshgrid(x, x). waterfall(x, y, z) Same as mesh but draws mesh lines in one direction only. 5 -95

Plots of the surface z xe (x 2 y 2) created with the mesh function and its variant forms: meshc, meshz, and waterfall. a) mesh, b) meshc, c) meshz, d) waterfall. Figure 5. 8– 4 5 -96