% urn.m
% This MATLAB code simulates an urn with two-types
% of balls (1 and 2). You are prompted to supply
% the initial number in the urn, the number of balls
% at which to stop simulating a particular urn, and
% the total number of urns to simulate. The rules are:
% choose a ball, replace the ball with itself and one
% other of the same type, repeat process.
%
% The output will be vectors: x contains a list of the
% final fraction of type1 balls in each urn; xct1 contains
% a list of the final fraction of type 1 balls in each urn for
% which the first ball chosen was of type1; and xct2 contains
% a list of the final fraction of type 1 balls in each urn for
% which the first ball chosen was of type2.
%
!c:
n1=input('how many number1 balls to start:    ')
n2=input('how many number2 balls to start:     ')
n=input('how many balls to stop urn at:     ')
q=input('how many simulations to run:     ')
kj=0;
kjj=0;
for i=1:q
ball1=n1;
ball2=n2;
en=n-n1-n2;
a=rand(1,1)
if a <= ball1/(ball1+ball2)
   ball1=ball1+1;
   ct=1;
else
   ball2=ball2+1;
   ct=2;
end;
for j=2:en
a=rand(1,1);
if a <= ball1/(ball1+ball2)
   ball1=ball1+1;
else
   ball2=ball2+1;
end;
end
num1(i)=ball1;
x(i)=ball1/(ball1+ball2);
if ct == 1
   kj=kj+1;
   xct1(kj)=x(i);
else   
   kjj=kjj+1;
   xct2(kjj)=x(i);
end
disp('number of ball1 and fraction is')
num1(i),x(i)
end
hist(x)
title('histogram for all urns'),pause
hist(xct1)
title('histogram for ball1 first urns'),pause
hist(xct2)
title('histogram for ball2 first urns'),pause