Team Lab 6
Matlab Programming: Loops and Logic
CS 310
In this team lab you will use Matlab programming to study some simple simulations. You will use loops and if statements to simulate a few random processes. This will give you practice in using loops and logical statements in programming.
Begin by starting Matlab on the computer.
1. Flipping Coins and Tossing Dice.
For any positive integer n the Matlab statement
randperm(n)
generates a sequence of the first n integers in a random order. (randperm is short for random permutation.)
1.a. Enter the command randperm(n) for several different values of n and see the results.
Use randperm(6) for tossing a die. The first value is the top of the die.
Use randperm(2) for tossing a coin. The first value of 1 is heads, otherwise tails.
1.b. The following Matlab loop simulates 20 tosses of a die and, for each toss, prints out the result.
tosses = 20;% This is how many times we will toss the die.
for toss= 1:tosses% This loops through all the tosses
value = randperm(6);% Get the value of the toss.
value(1)% Print out the top value.
end
A version that gives a nicer printout is the following:
clear ;
tosses = 20 ;% This is how many times we will toss the die.
% The disp function displays an array.
% In this case the array consists of a number (tosses) that is
% converted to a string, and a string (‘ tosses‘)
disp([ num2str(tosses), ' tosses'])
for toss = 1:tosses% This loops through all the tosses
value = randperm(6);% Get the value of the toss.
% disp displays an array with strings.
% Some of the strings are numbers that are converted to strings.
disp([' toss ',num2str(toss,2), ' value = ', num2str(value(1))])
end
Copy and paste either version of the code into a Matlab program file and run it several times.
Do you understand how the loop works?
Important!! Notice that we have chosen variable names that convey some meaning. We could have used ‘a’ instead of ‘toss’ and ‘b’ instead of ‘value’, and the code would have worked just as well. However, if you can not choose a name for a variable that conveys some sense of what that variable is used for, then you don’t know what you are doing!
1.c. We now want to roll the die 1000 times and see how many times we get a 3 on the die.
Copy and paste this code into Matlab and run it several times.
tosses = 1000;
count3 = 0 ;% This counts how many 3s we will get.
for toss = 1:tosses
value = randperm(6);
if value(1) == 3 % Check, is this toss is a 3?
count3 = count3 + 1 ;% If it is a 3, count it.
end
end
count3% Print out how many 3s we found.
Notice that == is used for ‘is the same as’ and = is used for assigning a value to a variable.
Notice that we need two end statements, one for the if statement and one for the for loop.
How many times did you get a 3 on the die?
What is the expected number of 3s in 1000 tosses?
Does the simulation agree with your intuition?
1.d.Nested Loops.
Run many trials of the simulation by putting the previous code inside another loop.
You may want to make a few small changes to make the code run nicely.
Trials = 100 ;% This is how many trial runs we will do.
Totals = 0 ;% This adds up the values of count3.
for trial = 1:Trials
( Put previous code here. )
Totals = Totals + count3 ;
end
Avg3 = Totals/Trials % This is the average number of 3s per trial.
This code is an example of nested loops, one inside the other.
Use ‘Smart Indent’ under the ‘Text’ menu to make your Matlab code look nicer.
If you do 100 trials what is the average number of 3s in 1000 tosses?
Is the computed result reasonable?
1.e. Breaking out of a loop.
Often we want to stop a loop when something special happens. For example, consider finding how many tosses of a die it takes on average to get ten 3s. In this case we can run the loop and when we get ten 3s we break out of the loop and record the number of the toss on which we got the tenth 3.
% Toss a die to determine how many tosses it takes to get ten 3's.
clear
% Set the parameters and initialize the total.
Trials = 1000;
tosses = 1000;
Total = 0 ;% For getting the average tosses to get 10 3s.
for trial =1:Trials
count3 = 0 ; % For each trial, reset the counter.
for toss = 1: tosses % For each trial do a bunch of tosses.
value = randperm(6) ;
if value(1) == 3
count3 = count3 + 1 ;
if count3 == 10 % Check to see if there are 10 3s.
break % If so, break out of the loop.
end
end
end
% The control goes here after the break.
Total = Total + toss ;% Add in the time to get 10 3s.
end
AvgTen3s = Total/Trials
What is the expected number of tosses that it takes to get ten 3s?
Do the computations agree with your intuition?
2. Random Walks.
2.a. Consider the problem of a bug hopping on a grid of equally spaced points. At any point on the grid, the bug has equal probability of going to the point to the right or the left. The next piece of code has the bug walking for a number of steps (given by tosses) and then printing out where the bug is at the end. This also shows the use of the else as part of the if statement. The plot and pause statements help show the position of the bug.
Run this code several times and see how it works.
% Random walk of a bug on a grid.
clear% Clear all variables.
clf % Clear the figure.
% Set the initial parameters.
tosses = 50;% Number of tosses.
Length = 10 ; % For how many points to plot.
position = 0 ; % Start in the middle
% Each toss is one time interval, so we now use the name time.
for time = 1: tosses
value = randperm(2) ;% Toss to see whether to go + or -.
if value(1) == 1
position = position + 1 ;% Right
else
position = position - 1 ;% Left
end
% Plot the current position and the grid.
plot( position, 0 ,'o', -Length:Length, 0 ,'+')
axis([ -(Length+1), Length+1, -.2 , .2])
% Put some text on the plot to show when and where.
text( -2, .1, ['time = ' num2str(time)] );
text( -2,-.1, ['position = ' num2str(position)]);
pause(.3)% Pause long enough to show what is happening.
end
position
Do you understand how this code works?
2.b. A Finite Grid. Falling off the end.
Consider now the case with only a finite number of grid points, indexed by the integers from –5 to 5. You are to find the average amount of time it takes for the bug to fall off the end of the grid if it starts in the middle. That is, if the bug gets to a grid point other than those in [-5,5], the trial is over and you have to start another trial.
Make the following changes to the previous code:
Reset Length to be 5.
Reset tosses to be 500.
Put in a break statement to stop the computation when the bug falls off the grid.
The break statement may look something like this:
if abs(position) > Length
break
end
Where do you put these statements?
Run the modified code several times.
2.c. Nest the previous code in another loop for several trials to get an estimate for the number of steps it takes before the bug falls off.
(If the bug starts at the center the expected number of time steps before falling off is (Length+1)2.)
3. Two-Dimensional Random Walks
Here we consider a bug on a two dimensional grid as shown. At each grid point, the bug has an equal chance of going in each of the four directions. The next piece of code generates a random walk for a bug starting at (0,0) on a grid. It also displays in a figure, the path and location of the bug.
3.a. Run this code several times and look at the code to see how it works.
(There is a figure here. Use the proper viewing option to see it.)
%Random walk of a bug in two dimensions.
clear
clf
MaxTime = 50;% Number of steps for the bug to move.
% This stuff is to plot a circle of radius Radius (for part 3.b.).
Radius = 5 ;% This the radius of the circle.
theta = 2*pi*(0:.02:1);% These are the angles to plot.
xcirc = Radius*cos(theta);% The x coordinates of the circle.
ycirc = Radius*sin(theta);% The y coordinates of the circle.
xpos = 0;% Initial x coordinate
ypos = 0;% Initial y coordinate
% Plot the initial position of the bug.
plot(xcirc, ycirc, xpos, ypos,'o' );
axis equal% Equal units in both directions.
pause(0.3)% Wait for a bit.
% xhist and yhist are for plotting the position history of the bug.
xhist(1) = xpos ;
yhist(1) = ypos ;
for time = 1: MaxTime
value = randperm(4) ;% Toss to find which way to go.
if value(1) == 1 % East
xpos = xpos + 1 ;
elseif value(1) == 2 % West
xpos = xpos - 1 ;
elseif value(1) == 3 % North
ypos = ypos + 1 ;
else % South
ypos = ypos - 1 ;
end
xhist(time+1) = xpos ;% Record the position.
yhist(time+1) = ypos ;
% Plot the history, the circle, and current position.
plot( xhist, yhist, xcirc, ycirc, xpos, ypos,'o' );
axis equal% Equal units in both directions.
% The text statement puts text at the specified point on the
% plot.
text( -6, 4, ['time = ' num2str(time)] );
text( -6, -4, ['position = ' num2str([xpos, ypos])]);
pause(0.3) % pause for a little while.
end
3.b. Can you determine the average time it takes for the bug to fall off the edge of the circle?
How you will check to see if the bug falls out of the circle?
If you plot inside of the loop for the trials you may need the statement
clear xhist yhist
to clear the history of the positions before starting the next path. You may want to ‘comment out’ all of the plotting to speed things up.
Extra Problems
X1. How would change the code for the one-dimensional bug walk if the probabilities were 1/3 to the right, 1/3 to the left, and 1/3 to stay put?
X2. How would change the code for the two-dimensional bug walk if the probabilities were 1/5 for each direction and 1/5 to stay put?