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

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)).