Skip to main content

Section 29.3 Beyond curve-fitting: source location

Suppose someone took water samples at various points \((x_k, y_k)\) and measured the concentration of some pollutant (microplastics, oil, etc) at each point. Let's use \(z_k\) for these amounts. A question of interest is the source of this pollution.

The theory of partial differential equations tells us that if a substance diffuses from a point source at \((u, v)\) under somewhat idealized conditions (ignoring currents, decay, etc) then its concentration after some time will be expressed by the two-dimensional Gaussian function

\begin{equation} z = A\exp\left(-((x-u)^2 + (y-v)^2)/s\right)\label{eq-2d-gaussian}\tag{29.3.1} \end{equation}

Here \(A, s\) are the shape parameters of the Gaussian (describing how tall and wide it is) and \((u, v)\) are its location parameters (where it is centered). We are really interested in the latter, but all four will need to be determined by the nonlinear least squares process.

One can visualize the above as surface-fitting: given the points \((x_k, y_k, z_k)\) we seek to fit a surface of type (29.3.1) to them, with the ultimate goal of locating the center of that Gaussian. Suppose the data is as follows.

x = [1 3 5 1 2 4 6 1 2 5 7];
y = [1 2 1 3 2 4 3 4 4 5 3];
z = [3 9 8 4 6 9 8 3 6 6 5];

Fit a two-dimensional Gaussian to the above data and report the expected location of the source.

Solution

The problem is solved in three lines: model function, optimization, output. The components of vector b below represent the parameters \(A, C, u, v\) in this order.

f = @(x, y, b) b(1)*exp(-((x-b(3)).^2 + (y-b(4)).^2)/b(2));
opt = fminsearch(@(b) sum((z - f(x, y, b)).^2), [1; 1; 1; 1]);
fprintf('The source is located near (%g, %g) \n', opt(3:4));