Page 4

Matlab: Buffon's Needle Problem

throws = 10000;   % percent makes the rest of the line a comment
x=rand(1,throws); % a vector of 10000 pseudorandom numbers
                  % in the range [0,1)
theta=rand(1,throws); % semicolon supresses printing
theta=0.5*pi*theta; % theta now 10000 random numbers between 0 and pi/2
hits= x<=0.5*sin(theta); % hits is a vector of 0's where false
                             % and 1's where true
sum(hits) % no semicolon prints the number of hits out of the 10000 throws
piEst=throws/sum(hits)

Page 6

Matlab: Histogramming Buffon's Needle Problem

throws = 100; % 100 trials
replications=2000;
for i=1:replications
  x=rand(1,throws);
  theta=0.5*pi*rand(1,throws);
  hits= x<=0.5*sin(theta);
  y(i)=sum(hits)/throws; % y= result vector
end
hist(y) % let hist pick the bins
recipEst=sum(y)/replications
vVec=(y-recipEst).*(y-recipEst); % vector of sqr deviations
v=sum(vVec)/(replications-1); % sample variance
stddev=sqrt(v) % sample standard deviation
piEst=1/recipEst

Page 9

Matlab: Gambler's Ruin

R=zeros(1,2000000); % vector of 2,000,000 zeros
i=1; R(i)=100; % gambler's initial fortune
while( R(i) > 0 & R(i) < 2100 )
 i=i+1;
 W = (rand  <  0.5); % random value of 0 or 1
 W = 2*W - 1; % random value of +/-1
 R(i)= R(i-1)+W; % gamblers new fortune
end
plot(1:length(R),R) % plot R against its index
 %%% run for investigating running time
for j=1:20
i=1; R(i)=100;
while( R(i) > 0 & R(i) < 2100 )
 i=i+1;
 W = 2*(rand  <  0.5)-1;
 R(i)= R(i-1)+W;
end
T(j)=i
end
hist(T)

Page 13

Matlab: Dice Rolls

red = floor(6*rand(1,1000))+1; % vector of 1000 random values 1,2,...,6
green  = floor(6*rand(1,1000))+1; % ditto for 2nd die
dice = red+green;
x=0:0.1:13; % vector of values from 0 to 13 incrementing by 0.1
for i=1:length(x)
 cdf(i)=sum(dice <= x(i));
end
plot(x,cdf)

Page 31

Matlab: Middle Sine Random Number Generator

tenpow5 = 100000;
x=0.5;
for i=1:10
  x=tenpow5*sin(x);
  w = floor(x); % first 5 digits
  x = x-w % digits after the 5th
end
 % now do some timing
tic
for i=1:100000
  x=tenpow5*sin(x);
  w = floor(x);
  x = x-w;
end
toc %0.5314
 % compare with Matlab's built-in generator
tic;
x=rand(1,100000);
toc; %0.0193

Page 32

Matlab: Middle Square Random Number Generator

tenpow7 = 1000000;
x = 0.1234567; % initial seed
for i=1:10
  x = x*tenpow7; % x now an integer
  w = x*x;  % square
  w = floor(w/1000); % drop last three digits
  w = w/tenpow7;  % 7 decimal digits
  x = w - floor(w) % drop the integer digits

Page 39

Matlab: Pi Estimation

nTrials = 10000;
x=rand(1,nTrials);
y=rand(1,nTrials);
hits=sum(x.^2+y.^2 < 1)
piEst = 4*hits/nTrials

Page 40

Matlab: Coupon Collecting

nLetters = 9; %ANGELFISH
nTrials = 10000;
for i=1:nTrials
    success = 0;
    nTries(i) = 0;
    for j=1:nLetters
        ANGELFISH(j)=0; %reset to letter not achieved
    end
    while success == 0
        nTries(i) = nTries(i)+1; %inc. count
        buy = 1+floor(nLetters*rand); %letter obtained
        ANGELFISH(buy) = 1;
        if sum(ANGELFISH)==nLetters
            success = 1;
        end
    end
end
hist(nTries)

Page 42

Matlab: Card Shuffling

n = 4; % a permutation of 4 objects
m = n; % do the whole deck
perm = 1:n %identity permutation
for i=1:m % if m=n the last is superfluous
    j = i+floor((1+n-i)*rand); % from the remaining
    swap = perm(i);
    perm(i) = perm(j);
    perm(j) = swap;
end
perm

Page 47

Matlab: Exercise 27.

low=0; high=0;
for i=1:10000000
  x=rand;
  if x < 0.1
    low = low + 1;
    y=rand; % for baseline comment this out
  else
    high = high + 1;
    y=rand; z=rand; % for baseline comment out
  end
end
low
high