% Control3_1.m
% Solution to problem 1 of control 3
% Jose M. Soler, May 2019

% Set parameters
clear
n = 100;  % total number of balls
b0 = 50;  % initial number of black balls
nx = 10;  % number of extractions per experiment
ne = 1e5; % number of experiments

% Perform extraction's experiment
c = zeros(nx,ne);                  % color of each extraction
for ie = 1:ne                      % loop on experiments
  cc = zeros(1,n);                 % inital ball colors
  cc(1:b0) = 1;                    % [0,1] <=> [white,black]
  for ix = 1:nx                    % loop on extractions
    k = randi(n-ix+1);             % index of extracted ball
    c(ix,ie) = cc(k);              % store color of extracted ball
    cc = [cc(1:k-1),cc(k+1:end)];  % exclude extracted ball
  end
end

% Find and plot probability of extracting b black balls
be = sum(c);            % number of black balls in each experiment
for b = 0:nx            % loop on possible number of extracted black balls
  p(b+1) = numel(find(be==b))/ne;
end
figure(1)
bar(0:nx,p)
xlabel('Number of black balls extracted')
ylabel('Probability')

% Probability of seventh extracted ball being white (not black)
y = ~c(7,:);                         % is 7th ball white?
ny = numel(y);                       % number of values
y0 = sum(y)/ny;                      % mean value (probability of true)
dy = sqrt(sum((y-y0).^2)/ny/(ny-1)); % standard dev. of mean value
fprintf('Probability of 7th ball being white = %6.4f +/- %6.4f\n',y0,dy)

% Probability of two succesive black balls
y = c(1:nx-1,:).*c(2:nx,:);          % were succesive balls black?
ny = numel(y);
y0 = sum(y(:))/ny;
dy = sqrt(sum((y(:)-y0).^2)/ny/(ny-1));
fprintf(' Prob. of two succesive black balls = %6.4f +/- %6.4f\n',y0,dy)
