1 / 9

Physics 114: Lecture 19 Least Squares Fit to 2D Data

Physics 114: Lecture 19 Least Squares Fit to 2D Data. Dale E. Gary NJIT Physics Department. Nonlinear Least Squares Fitting. We saw that for a function y ( x ) the chi-square is which can be used to fit a curve to data points y i .

jstoudt
Télécharger la présentation

Physics 114: Lecture 19 Least Squares Fit to 2D Data

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Physics 114: Lecture 19 Least Squares Fit to 2D Data Dale E. Gary NJIT Physics Department

  2. Nonlinear Least Squares Fitting • We saw that for a function y(x) the chi-square is which can be used to fit a curve to data points yi. • In just the same way, we can easily extend this idea to fitting a function z(x,y), that is, a surface, to a 2D data set of measurements zi(x,y). • You know that it takes three parameters to fit a Gaussian curve • To extend this to a 2D Gaussian is straightforward. It is: • Note that this has 5 parameters, but the most general way to describe a 2D Gaussian actually has 6 parameters, which you will do for homework.

  3. Plotting a 2D Gaussian in MatLAB • Let’s introduce some useful MatLAB commands for creating a 2D Gaussian. • One can use a brute force method such as • z = zeros(101,101); • a = 2.5; bx = 2.3; by = -1.2; cx = 1.5; cy = 0.8; • i = 0; • for x = -5:0.1:5 • i = i + 1; • j = 0; • for y = -5:0.1:5 • j = j + 1; • z(i,j) = a*exp(-((x-bx)/cx)^2 – ((y-by)/cy)^2); • end • end • imagesc(-5:0.1:5,-5:0.1:5,z) • Or, we can plot it as a surface to provide x and y coordinates • surf(-5:0.1:5,-5:0.1:5,z)

  4. Better Way of Creating 2D Functions • That way of using for loops is a very slow, so here is another way of doing it for “round” Gaussians: • x = -5:0.1:5; • y = -5:0.1:5; • siz = size(x); siz = siz(2); • a = 2.5; bx = 2.3; by = -1.2; c = 0.8; • dist = zeros(siz,siz); • for i=1:siz • for j=1:siz • dist(i,j) = x(i)^2 + y(j)^2; % Create “distance” array once • end • end • % use the next two statements multiple times for different Gaussians • sdist = circshift(dist,round([bx/0.1,by/0.1])); % Shifted “distance” array • z = a*exp(-sdist/c); • imagesc(x,y,z);

  5. Best Way of Creating 2D Functions • The best way, however, is to use the meshgrid() MatLAB function: • [x, y] = meshgrid(-5:.1:5, -5:.1:5); • a = 2.5; bx = 2.3; by = -1.2; cx = 1.5; cy = 0.8; • z = a*exp( - (((x-bx)/cx).^2 + ((y-by)/cy).^2) ; • surf(x,y,z); • This is completely general, and allows vector calculations which are very fast. • Here is another example, the sinc() function: • z = a*sinc(sqrt(((x-bx)/cx).^2 + ((y-by)/cy).^2)); • You can see that the meshgrid() function is very powerful as a shortcut for 2D function calculation. • To add noise, just add sig*randn(size(z)), where sig is the standard deviation of the noise.

  6. Searching Parameter Space in 2D • Just as we did in the 1D case, searching parameter space is just a trial and error exercise. • % Initial values around which to search • a0 = max(max(z)); • [ix iy] = find(z == a0); • bx0 = x(ix,iy); by0 = y(ix,iy); cx0 = 1; cy0 = 1; • astep = 0.05; bstep = 0.05; cstep = 0.05; • chisqmin = 1e20; • for a = a0-5*astep:astep:a0+5*astep • for bx = bx0-5*bstep:bstep:bx0+5*bstep • for by = by0-5*bstep:bstep:by0+5*bstep • for cx = cx0-5*cstep:cstep:cx0+5*cstep • for cy = cy0-5*cstep:cstep:cy0+5*cstep • ztest = a*exp(-(((x-bx)/cx).^2 + ((y-by)/cy).^2)); • chisq = sum(sum(((ztest-z)/sig).^2)); • if (chisq < chisqmin); chisqmin = chisq; params = [a,bx,by,cx,cy]; end • end • end • end • end • end params = 2.1465 2.3000 -1.2000 1.2500 0.8000

  7. The Result Looks Good, but… • To calculate the result of the fit, just type • ztest = params(1)*exp(-(((x-params(2))/params(4)).^2 + ((y-params(3))/params(5)).^2)); • imagesc(x(1,:),y(:,1),ztest) • This gives a reasonable result, but compare the fit parameters with the (known) params: • a = 2.1; bx = 2.3; by = -1.2; cx = 1.3; cy = 0.8; • params = 2.15 2.3 -1.2 1.25 0.8 • The value of cx is especially far off, because our search window was limited to 0.75-1.25, so it could not reach 1.3. • Look at the residuals: • imagesc(x(1,:),y(:,1),ztest-z) • Can you see a pattern in the noise? • To bring out the pattern better, let’s introduce convolution.

  8. 1D Convolution and Smoothing • Let’s create a noisy sine wave: • u = -5:.1:5; • w = sin(u*pi)+0.5*randn(size(u)); • plot(u,w) • We can now smooth the data by convolving it with a vector [1,1,1], which does a 3-point running sum. • wsm = conv(w,[1,1,1]); • whoswsm Name Size Bytes Class Attributes wsm 1x103 824 double • Notice wsm is now of length 103. That means we cannot plot(u,wsm), but we can plot(u,wsm(2:102)). Now we see another problem. • Try this: • plot(u,w) • hold on; plot(u,wsm(2:102)/3,’r’)

  9. 2D Convolution • To do 2D smoothing, we will use a 2D kernel k = ones(3,3), and use the conv2() function. So to smooth the residuals of our fit, we can use • zsm = conv2(ztest-z,k)/9.; • imagesc(x,y,zsm(2:102,2:102)) • Now we can see the effect of missing the value of cx by 0.05 due to our limited search range. • There are other uses for convolution, such as edge detection. For example, we can convolve with a kernel k = [1,1,1,1,1,1]. • Or a kernel k = [1,1,1,1,1,1]’. Or even a kernel k = eye(6). Or k = rot90(eye(6)).

More Related