% Control3_1.m
% Solution to problem 1 of control 3. 2023/05/05
% Jose M. Soler

clear
nb = 5;       % number of extracted balls
niter = 1e4;  % number of iterations (experiments)

% Set ball values
xb0 = [ 1*ones(1,16), 2*ones(1,13), 3*ones(1,11) ];

% Loop on experiments
for iter = 1:niter
  xb = xb0;                           % each experiment begins with all balls
  for ib = 1:nb                       % loop on extractions
    jb = randi(numel(xb));            % index of extracted ball
    x(ib,iter) = xb(jb);              % value of extracted ball
    xb = [xb(1:jb-1),xb(jb+1:end)];   % remove extracted ball from list
  end
end

% Alternatively
% load('BolasExtraccion.mat');
% x = BolasExtraccion;

% Plot histogram of sums of extracted balls
figure(11)
histogram(sum(x),'Normalization','Probability');
xlabel('Sum of extracted ball values')
ylabel('Probability')
grid on
set(gca,'FontSize',15)

% Print some probabilities
fprintf('               P(x(4)==3) = %5.3f\n', sum(x(4,:)==3)/niter )
fprintf('            P(~any(x==1)) = %5.3f\n', sum(~any(x==1))/niter )
fprintf('P(any(x(i:i+2)==[1,2,3])) = %5.3f\n', sum( all(x(1:3,:)==[1;2;3]) | ...
                                                    all(x(2:4,:)==[1;2;3]) | ...
                                                    all(x(3:5,:)==[1;2;3]) )/niter )
