Suppose it's night time, and you need to estimate the area of a lake that you can't see. All you know is that it lies inside a square park that is a mile on each side. Oh, and you have a cannon. One that is so inaccurate that if you fire it into the park, the cannonball will land somewhere in the park completely at random. How do you estimate the area of the lake?
Well here's a plan. Let's say you fire once. You don't hear a splash, so you estimate the area of the lake as 0 square miles. You fire again, and you hear a splash. So now you estimate that the lake is 1/2 square mile. After firing 100 times, and hearing 47 splashes, you guess it's about 0.47 square miles. If you were very patient, it's reasonable to assume you could actually make a very good estimate of the area of the lake.
Suppose your park is the square
-5 <= X <= 5,
-5 <= Y <= 5;
and that the lake is the ellipse
9 x2 + 16 y2 = 144.
Draw a picture of the region, and estimate the area by sight.
MATLAB does not have a good facility for plotting an implicit relationship like this. If we really want a nice picture of the ellipse (or any other implicit function) we have to be a little sneaky. First, define a function:
Z(X,Y) = 9 x2 + 16 y2 - 144.
This function is negative inside the ellipse, zero on the ellipse,
and positive outside. We need to define a set of point (X,Y)
that are evenly spaced in our range, and then evaluate Z there.
Similar to linspace for a 1 dimensional range, we now use
[x y] = meshgrid ( -5:0.2:5, -5:0.2:5 )
to generate a grid of points. Now we just have to issue a command
that defines the value of Z.
Once we have the Z array set up, MATLAB can draw a contour plot with this command:
contour ( z )
Unfortunately, we see much more than we wanted: lots and lots of
ellipses. We can't even be sure the ellipse we want is one of them.
Luckily, we can tell MATLAB which ellipse to draw by specifying the
contour levels. Because we only want one contour level, the command
looks a little strange:
contour ( z, [0 0] )
If this doesn't quite work, try expanding the range of values that
will show up:
contour ( z, [-0.001 +0.001] )
The ellipse will not be its true shape, because MATLAB's graph paper is not square. To fix this, type:
axis square
Also, MATLAB doesn't know what the X and Y coordinates are. If we
want them to show up on the axes, we need to add them in the command:
contour ( x, y, z, [0 0] )
If we use the grid on command, we'll even get grid lines
that help us to estimate the area.
Exercise A: from the hints above, create a plot of the ellipse. Use the grid to make a crude guess for the area of the ellipse. When you give your report, the line drawing of the ellipse won't show up well. Try the contourf command, which does the same job, but "fills in" the region.
Now we're ready to try to estimate the area of the lake using simulated cannon fire. We'll fire at the park N times, each time landing at a random point (X,Y) where we will evaluate Z. Since any point inside the lake will have a negative Z, we can estimate the area of the lake by multiplying the total area of the park by the fraction of the Z values that are negative.
Our MATLAB M file might look like:
function area = cannon ( n )
but you might at first be puzzled about how to write down the
steps we need to take.
Here are two important tools that can help you make this computation. For the array A:
1 2 3
4 5 6
we can make an array of the same shape, containing 0's and 1's,
corresponding to where the array is less than 5, as follows:
B = ( A < 5 )
which results in
0 0 0
0 1 1
The second tool is the SUM function. If X is a vector, then SUM(X) is the sum of the elements. If A is an array, then SUM(A) is a row vector containing the sum of the elements of each column.
Believe it or not, these clues should be enough to allow you to make your entire are computation in about 5 commands.
Exercise B: fill in the details for the routine cannon.m. Compute area estimates using 1, 10, 100, 1000, 10,000 and 100,000. points. Make a table of the number of points, the area estimate, and the error. Is the absolute value of the error going down? Can you estimate how fast it is going down? Is there a way to check whether your estimate is good?
Now it would be nice to see a plot of the points that we are using. But how are we going to do that? Here's an idea...The MATLAB PLOT command can be told to just plot points, no lines between them, if we follow the X and Y coordinates by a style argument that has the letter "o":
plot ( x, y, 'o')
so an idea here would be to draw the dots as we find them.
A first draft for a routine to do this might look like this:
function cannon_plot ( n )
x = ???
y = ???
z = ???
plot ( x, y, 'o' )
Try this. Use N = 10, 100, 1000. It's sort of a mess.
This is half of a good idea, but it's not really giving us the information we need. What what be much better would be to see which dots hit the lake and which don't. I see blue lake dots, and green park dots. How can we do this? Well, the color is easy. Typing
plot ( x, y, 'bo')
will cause the dots that are plotted to be blue. But what I need
is a command that says:
plot ( x(LAKE), y(LAKE), 'bo', x(PARK), y(PARK), 'go' )
where LAKE is somehow telling MATLAB to only pick the points that
go into the lake.
Luckily, there's another MATLAB command that will do that for us. The FIND command will make a list of all the entries of a vector or array for which something is true. We can then use that list to select just those elements. For instance:
x = [ 2, 7, 4, 0, 9, 3, -10, 6 ]
list = find ( x < 5 )
x(list)
After the second command, list contains the indices 1, 3, 4, 6 and
7. And x(list) will be exactly the values [2,4,0,3,-10].
Exercise C: Using the idea of a list, and the information in Z, revise the cannon_plot command so that it plots blue lake dots and green park dots. If you get this right, then a plot for N = 1000 should actually look nice.
What we were calling the cannonball idea is really called random sampling. We try to find out an overall property of a function or a region or an object by finding values of that property at points chosen at random. You might think there are better ways to measure the area of an ellipse than firing a cannonball at it, but suppose I had given you a very complicated shape with all sorts of crazy angles? Sometimes the method of random sampling is the only effective way of handling difficult problems.
Now we're going to take this idea of random sampling of a region, and use it for a more interesting problem, related to the idea of a Voronoi diagram. The Voronoi diagram is a very useful tool for building up an understanding of a whole area based on information at selected points. You might use a Voronoi diagram, for instance, if you had rainfall measurements at 10 spots throughout the county, and wanted to estimate how much rain fell overall, and how much rain fell at particular points where no data was collected.
Suppose that you have a garden that is a square 100 feet on a side, and that there are 10 very hungy ants, each of whom has a nest somewhere in the garden. From time to time, you throw a sugar cube into the garden, and all the ants run out of their nests towards the sugar cube. The ants run at the same speed, so whichever ant's nest is closest to the sugar cube is the one who will get it. This means that each ant has a territory, inside of which it is guaranteed to catch every sugar cube that lands there. The ants want to find out the shape of these regions, so the losers can stay in their nests and stop wasting time. Can we tell them?
The ant nests are at these points:
# X Y Area?
-- -- -- ------
1 100 80 ______
2 0 0 ______
3 25 65 ______
4 40 20 ______
5 94 20 ______
6 92 75 ______
7 30 20 ______
8 90 20 ______
9 45 2 ______
10 35 100 ______
We'll use the random sampling idea to compute the areas. Tossing in the sugar cube is just going to be selecting two random values between 0 and 100. The tricky part is that now we have to figure out which ant is closest to the point. We can make a variable called hits which keeps track of how many times a sugar cube landed in each ant's territory. Then, for ant number 3, for instance, the area of its territory is approximated by (hits(3)/n) * 10000.
Here is a sketch of the MATLAB M file ants.m that you will need.
function area = ants ( n )
nest_x = [ 100; 0; 25; ... ];
nest_y = [ 80; 0; 65; ... ];
hits = zeros(10,1);
x = rand ( n, 1 );
y = rand ( n, 1 );
for i = 1 : n
nearest = 1
dist_min = sqrt ( ( x(i) - nest_x(1) )^2 ...
+ ( y(i) - nest_y(1) )^2 );
for j = 2 : 10
dist = sqrt ( ( x(i) - nest_x(j) )^2 ...
+ ( y(i) - nest_y(j) )^2 );
if ( dist < dist_min )
nearest = j;
end
end
hits(nearest) = hits(nearest) + 1;
end
area(1:10) = 100.0 * 100.0 * hits(1:10) / n;
Exercise D: Get the ants.m program running. Try to find a large enough value of N that your results seem "accurate". About how many digits of your answers do you have confidence in?
Exercise E: To get an idea of the shape of the territories, we can try the same idea we used for the lake. This time, however, we have 10 sets of data to plot. I can see two ways of doing this:
i1 = find ( nearest == 1 );
i2 = find ( nearest == 2 );
...
i10 = find ( nearest == 10 );
plot ( x(i1), y(i1), 'ro', x(i2), y(i2), 'go', ...
x(i3), y(i3), 'bo', ...
Exercise F: For the ant's nest data, try the MATLAB command VORONOI, and compare it to the picture you computed.
Back to the Mini Projects page.
Last revised on 08 June 2001.