Team 2 - Segment 21 - Class Activity

From Computational Statistics Course Wiki
Jump to navigation Jump to search

Team 2:

Todd Swinson,

Elad Liebman,

Shuqi Zhang


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]);


Seg21 plot1.png


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,:);