7 - Admin Items Lecture 6: Nonlinear Curve Fitting Oct. 28...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: Admin Items Lecture 6: Nonlinear Curve Fitting Oct. 28 – Last “regular” lecture. Assignment 7 will Oct. regular” be posted. No lecture on Nov. 4, but the final assignment will No be posted (Assignment 8)! November 11 – Review and Research lecture. November Nov 18 – Study break for final quiz. No lecture. Nov November 25, 1-2:30pm – Final Quiz (in class, 90 (in November 1min). Worth 15%. Prof. M. Dobbs, Physics 257 October 25, 2007 2 Review: Curve Fitting χ2 = ∑ i ( yi − f ( xi ; ai ) )2 Δyi 2 Review: Curve Fitting χ = ∑ ( y − f ( x ; a )) σ 2 i i 2 i i 2 We assume the data points are drawn from a We Gaussian distribution. The probability of randomly getting one of the data points The is: ⎛ −( y − f ( x ) ) ⎞ P= 1 2πσ 2 e⎝ ⎜ ⎜ 2 i i So the probability of getting all of the data points is: So P=∏ i i 1 2πσ i 2 e ⎛ − ( yi − f ( xi ) )2 ⎞ ⎜ ⎟ ⎜ ⎟ 2σ i 2 ⎝ ⎠ 2σ i 2 ⎟ ⎟ ⎠ Let’s take the logarithm, and we get: Let’ Recall that for dice, when we require two outcomes Recall (such as rolling two 3’s), we multiply the 3’ probabilities, P = 1/6 x 1/6 = 1/36. So the probability of getting all of the data points is: So P=∏ i ln( P) = ∑ i 2 ⎛1 − ( yi − f ( xi ; ai ) ) + ∑ ln⎜ 2 2 ⎜ 2σ i i ⎝ 2πσ i log likelihood = 1 2πσ i 2 e⎝ ⎛ − ( yi − f ( xi ) )2 ⎞ ⎜ ⎟ ⎜ ⎟ 2σ i 2 ⎠ Where we have dropped the terms that don’t depend on Where don’ the fit parameters (just a constant offset). 3 − ( y − f ( x ; a )) 1 ∑ i σ2i i 2i i ⎞ ⎟ ⎟ ⎠ 2 Prof. M. Dobbs, Physics 257 October 25, 2007 Prof. M. Dobbs, Physics 257 October 25, 2007 4 Review: Standard Error from χ2 Fits The best fit parameters correspond to the minimum The χ2min of the χ2 function. 2 χ2 = ∑ i Nonlinear Curve Fitting: Intro The method we used previously to find the optimum value of The fit-parameters is limited to fitting functions of the type: fit- ( yi − f ( xi ; ai )) σi 2 y ( x) = ∑ a j f j ( x) j =1 m i.e. functions that are linear in the ‘m’ parameters to be i.e. determined (the aj above), since in this case minimizing ChiChisquare simultaneously in all parameters (by taking partial derivatives and setting them to zero) results in a set of ‘m’ coupled equations in ‘m’ unknowns. This is a simple exercise in linear algebra… problem solved! This algebra… Can be used for all linear functions, linearizable functions or Can linear combination of functions (e.g. polynomials). See handout accompanying lecture notes on WebCT for a discussion of fitting to higher order polynomials. Prof. M. Dobbs, Physics 257 October 25, 2007 5 Prof. M. Dobbs, Physics 257 October 25, 2007 6 1 Best Fit Parameters for Least-squares fitting Leastof data to a quadratic function Application: Estimating Slopes of Curves Fit a high order polynomial (this is Fit phenomenological – there is no theory here!): y = −0.0000 x 5 + 0.0013x 4 − 0.0230 x 3 + 0.1371x 2 + 0.1793 x − 0.1646 dy = −0.0000 x 4 + 0.0052 x 3 − 0.0690 x 2 + 0.2742 x + 0.1793 dx y ( x) = a1 + a2 x + a3 x 2 slope = The derivative of this equation is the The slope. This allows us to calculate the slope m This at any point x. We can plug in x=1 to find the slope at We the first data point. >> p=polyfit(x,y,5) p= -0.0000 0.0013 0.0230 0.1371 0.1793 -0.1646 >> line = polyval(p,x) polyval(p,x) >> plot(x,y,'.',x,line,'-') plot(x,y,'.',x,line,'Prof. M. Dobbs, Physics 257 October 25, 2007 8 Prof. M. Dobbs, Physics 257 October 25, 2007 7 Try it in Matlab: Example shown in lecture Matlab: y 5 4 3 2 1 0 -1 -2 -3 -4 -5 5 4 3 2 1 y2 0 -1 -2 -3 y = 0.05 x3 – 0.5 x2 +0.5 x 0 1 2 3 4 5 -4 5 4 3 x 6 7 8 9 10 -5 0 1 2 3 4 5 10 x 6 7 8 9 10 dy/dx (for data and fit) y2 and yfit 2 1 0 -1 -2 -3 -4 -5 yfit = 0.0501 x3 – 05142 x2 +0.6169 x – 0.1165 8 6 4 2 0 -2 -4 -6 -8 0 1 2 3 4 5 Prof. M. Dobbs, Physics 257 October 25, 2007 9 x 6 7 8 9 10 0 1 2 3 4 5 Prof. M. Dobbs, Physics 257 October 25, 2007 x 6 7 8 9 10 10 Nonlinear Curve Fitting: Intro We need to take a different approach for fitting to We functions that are nonlinear in the fit parameters. Common examples include Lorentzian functions, multiCommon multiexponential decay curves (and many, many more…) more… Nonlinear Curve Fitting There are several approaches that are often used in order to There search for the minimum. They all require Selecting starting values for the parameters Selecting Some algorithm for making changes to the starting values in order to Some order search out the global minimum value of Chi-Square. Chi- For such problems we must treat Chi-Square as a For Chicontinuous function of the ‘m’ parameters (this function describes a “hypersurface” in an mhypersurface” mdimensional space) and search that parameter space for the minimum value. Prof. M. Dobbs, Physics 257 October 25, 2007 11 The most obvious approach is “Brute Force”… just map out The Force”… the surface over a regular ‘grid’ of trail values that covers the grid’ important region of parameter space. There is some ‘art’ in art’ all of this…. this… More sophisticated approaches apply a ‘gradient search’ to More search’ search out minima and are often necessary when the number of parameters is large. The most often used algorithm is called the “Levenberg-Marquardt method” Levenbergmethod” Prof. M. Dobbs, Physics 257 October 25, 2007 12 2 Numerical Curve Fitting Example Create a file chi2_line.m Create containing the χ2 function for a straight line. Plotting χ2 We clearly see the minimum in χ2, which corresponds to the best fit line. χ 2line(x,y,1.57,a,b) 30 14 13 12 4.9 20 19 18 17 χ2 = ∑ i ( yi − f ( xi ; ai ))2 Δyi 2 x = [1:1:15]; y=[6.5,12.3,16.2,21.9,26.5,32. 3, 38.2,44.8,49.0,54.1,55.6,60 .4,64.7,74.0,77.5]; a=5.019 ; b=2.13 ; >> chi2_line(x,y,1.57,a,b) ans = 12.9469 >> ezsurfc( @(a,b) chi2_line(x,y,1.57,a,b), [4.9,5.1,1.5,2.5] ) >> fplot(@(a) chi2_line(x,y,1.57,a,b),[4. 9,5.1]) >> fplot(@(b) chi2_line(x,y,1.57,a,b),[1, 3]) χ2 16 15 function [chi2] =chi2_line(x,y,ey,a,b) % This function returns the chi2 for a straight line y=ax+b for error ey chi2 =sum( (y-a*x-b).^2 )/ey^2 end 25 20 4.92 4.94 4.96 4.98 5 5.02 5.04 5.06 5.08 5.1 a 21 20 19 15 2.5 χ2 2 18 17 16 15 14 13 12 Let’s do something new, and plot Let’ the χ2 function itself. Notice the clear minimum at Notice a=5.019 ; b=2.13 b 1.5 4.9 4.95 a 5 5.05 5.1 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 b Prof. M. Dobbs, Physics 257 October 25, 2007 14 Prof. M. Dobbs, Physics 257 October 25, 2007 13 Finding the χ2 Minimum For simple functions like a polynomial, it is easy to find the For minimum by simply differentiating the χ2 function with respect to each parameter, setting each equal to zero and solving the set of linear equations that result. We did this explicitly using matrix mathematics for a degree 1 We polynomial (straight line). Finding the χ2 Minimum fminsearch() fminsearch() p = (slope, y-intercept) (slope, yStarting values for p(1) and p(2) Starting >> chi2_line(x,y,1.57,a,b) ans = 12.9469 >> result = fminsearch(@(p) chi2_line(x,y,1.57,p(1),p (2)), [1,1] ) For more complicated functions, it is often difficult or For impossible to find the minimum analytically (as we did for the straight line). Instead, we can find it numerically, by asking the computer to Instead, walk around the χ2 curve until it finds the global minimum. In compu-speak, this is called minimization, and matlab In compuprovides this functionality through the fminsearch() function. fminsearch() Prof. M. Dobbs, Physics 257 October 25, 2007 15 This search for the minimum This value returns our best fit values for the line; i.e. the parameter ‘a’ and ‘b’ that minimize Chi-square. ChiWhat about the uncertainties in What these values? How doe we estimate those using numerical techniques. result = 5.0193 2.1124 Prof. M. Dobbs, Physics 257 October 25, 2007 16 Variation of Chi-Square near the minimum ChiLast lecture we determined the natural logarithm of the likelihood fcn’ for a given data set ‘yi’ and fit fcn’ evaluated with particular set of parameters {a1, a2,…,an}: Variation of Chi-Square near the minimum Chi- ln( P) = ∑ i 2 ⎛1 − ( yi − f ( xi ) ) + ∑ ln⎜ 2 2 ⎜ 2σ i i ⎝ 2πσ i ⎞ ⎟ ⎟ ⎠ Chi-square varies as the square of the distance from a minimum in each parameter. An increase of 1 standard deviation (σaj) in the parameter from the optimum value, aj’, increases Chi-square by 1. In general, variations of n* σaj increase Chi-square by ‘n’. This provides an easy way of estimating the uncertainties in the best-fit parameters. Prof. M. Dobbs, Physics 257 October 25, 2007 17 Prof. M. Dobbs, Physics 257 October 25, 2007 18 3 Plotting Standard Errors Standard errors for a best fit are Standard normally plotted as contour intervals for 1,2 and 3σ. Make a contour plot of the χ2 Make function. function. χ 2line(x,y,1.57,a,b) 4 Plotting Standard Errors χ 2line(x,y,1.57,a,b) 4 3.5 3 a=5.019 ; b=2.13 ; cf = ezcontourf( @(a,b) chi2_line(x,y,1.57,a,b), [4.8,5.2,0,4] ) chimin = 1σ ≈ 68.3% Error Contour 2σ ≈ 95% Error Contour 2.5 2 1.5 1 0.5 3.5 3 2.5 2 1.5 1 0.5 0 4.8 chi2_line(x,y,1.57,a,b) clines = [chimin, chimin+1, chimin+2, chimin+3 ] set(cf,'LevelList',clines) 3σ ≈ 99.7% Error Contour b 4.85 4.9 4.95 5 a 5.05 5.1 5.15 5.2 b 0 4.8 Prof. M. Dobbs, Physics 257 October 25, 2007 19 4.85 4.9 4.95 5 a 5.05 5.1 5.15 5.2 20 Prof. M. Dobbs, Physics 257 October 25, 2007 Plotting Standard Errors χ 2line(x,y,1.57,a,b) 4 3.5 3 2.5 Numerical Methods Notice thus far that we have developed a completely Notice generic method for fitting any function and determining the error contours. χ2min Δb = 0.85 Best fit parameters We don’t really need to know anything about the function We don’ other than how to evaluate it. b = 2.13 2 1.5 1 This is a powerful technique, and we can apply it to This any function we care to dream up. Δa = 0.094 b 0.5 0 4.8 a = 5.019 4.85 4.9 4.95 5 a 5.05 5.1 5.15 5.2 21 Sometimes it is not so easy to determine the ‘global Sometimes minimum’.. minimum’ Prof. M. Dobbs, Physics 257 October 25, 2007 22 Prof. M. Dobbs, Physics 257 October 25, 2007 Summary Numerical minimizations of Numerical χ2 = ∑ i ( yi − f ( xi ; ai ))2 Δyi 2 yields the best fit parameters for the function f; i.e. the set {a1, a2,…,an} that minimize chi-square. a2,… chiIn Matlab, this can be done with the fminsearch() function. In fminsearch() The χ2 function can be plotted, and the confidence The interval contours superimposed. The 1σ uncertainty occurs at χ2 = χ2min + 1 The The nσ uncertainty occurs at χ2 = χ2min + n The Prof. M. Dobbs, Physics 257 October 25, 2007 Prof. 23 4 ...
View Full Document

Ask a homework question - tutors are online