Team 2 - Segment 21 - Class Activity
The background for this activity can be found here: Find the Volcano.
We are given a model with 5 unknown parameters that vary with each eruption.
Since we are interested in which area of the country to evacuate, we want to estimate . So, we create a chi-square function using the model, and attempt to minimize it. In other words, find the parameters that minimize the chi-square function.
Using what we can infer from the measurements, we provide an initial guess of the parameters to help guide the minimization algorithm. We believe the noisy measurements place the eruption hot spot around and we noticed that the max temperature was around 160 while the minimum temperature was near 0, so .
We now attempt to minimize our chi-square function with MATLAB's fminsearch routine.
ymodel = @(x, y, b) b(1) + (b(2) - b(1)) * exp(-1 * ((x - b(3)).^2 + (y - b(4)).^2) / ( 2 * b(5)^2 )); chisqfun = @(b) sum( ((t - ymodel(x, y, b)) ./ 30).^2); bguess = [0 160 60 60 30]; bfit=fminsearch(chisqfun, bguess) bfit = 17.6735 195.3020 56.0662 79.0983 38.8913
The estimate for is which seems reasonable given the measurement plots from the problem description page.
But now we must provide the authorities with a region to evacuate, as they can not evacuate a single point. In other words, what is the error ellipse of our location estimate?
We start by drawing an ellipse at 1 standard deviation around the eruption point. (We will assume that our estimates are drawn from a multivariate Normal distribution.)
To do this, we use the Hessian method to obtain the covariance matrix for the MVN of our estimate. MATLAB has a handy function that will produce a Hessian matrix corresponding to the minimization of a function:
[X,FVAL,EXITFLAG,OUTPUT,GRAD,HESSIAN] = fminunc(chisqfun, bguess); cov = inv(0.5*HESSIAN) cov = 1.0e+03 * 0.3823 -0.0353 0.0175 0.0699 -0.0833 -0.0353 1.6430 -0.0753 -0.3470 -0.2159 0.0175 -0.0753 0.0130 0.0204 0.0064 0.0699 -0.3470 0.0204 0.1013 0.0323 -0.0833 -0.2159 0.0064 0.0323 0.0519
Now that we have computed the covariance matrix, we need to marginalize out the parameters we are not interested in. Here, we keep only , which are the 3rd and 4th parameters in our parameter vector.
marg = cov(3:4, 3:4) marg = 13.0206 20.4205 20.4205 101.2641
Now we use the errorellipse function provided in Segment 17 to draw the error ellipse at 1 standard deviation around the eruption point.
mu = bfit(3:4)'; [a b] = errorellipse(mu, marg, 1, 100); plot(a, b); axis([0 100 0 100]);
We now have a reasonable estimate of where the eruption will take place, and can alert the authorities. However, it might be more useful to provide them with a 99% confidence interval (which will be left for the next class).
The full MATLAB code is provided below:
ymodel = @(x, y, b) b(1) + (b(2) - b(1)) * exp(-1 * ((x - b(3)).^2 + (y - b(4)).^2) / ( 2 * b(5)^2 )); chisqfun = @(b) sum( ((t - ymodel(x, y, b)) ./ 30).^2); bguess = [0 160 60 60 30]; bfit=fminsearch(chisqfun, bguess) [X,FVAL,EXITFLAG,OUTPUT,GRAD,HESSIAN] = fminunc(chisqfun, bguess); cov = inv(0.5*HESSIAN); marg = cov(3:4, 3:4); %marginalize out the parameters we aren't interested in mu = bfit(3:4)'; [a b] = errorellipse(mu, marg, 1, 100); plot(a, b); axis([0 100 0 100]);
function [x y] = errorellipse(mu,sigma,stdev,n) L = chol(sigma,'lower'); circle = [cos(2*pi*(0:n)/n); sin(2*pi*(0:n)/n)].*stdev; ellipse = L*circle + repmat(mu,[1,n+1]); x = ellipse(1,:); y = ellipse(2,:);