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, 12: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 fitparameters 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 Leastsquares 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 ChiSquare. Chi For such problems we must treat ChiSquare 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 “LevenbergMarquardt 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( (ya*xb).^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, yintercept) (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 compuspeak, 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 Chisquare. 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 ChiSquare 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 ChiSquare near the minimum Chi ln( P) = ∑
i 2 ⎛1 − ( yi − f ( xi ) ) + ∑ ln⎜ 2 2 ⎜ 2σ i i ⎝ 2πσ i ⎞ ⎟ ⎟ ⎠ Chisquare 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 Chisquare by 1. In general, variations of n* σaj increase Chisquare by ‘n’. This provides an easy way of estimating the uncertainties in the bestfit 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 chisquare. 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
 Fall '07
 Dobbs
 Regression Analysis, Prof. M. Dobbs

Click to edit the document details