% Sum of Bernoulli variables is binomial
% Textbook: Part I, "Discrete Random Variable", Bernoulli and binomial distributions
%
% X = sum of n independent Ber(p) draws, so X ~ Bin(n,p) with
% p_X[k] = nchoosek(n,k) p^k (1-p)^(n-k), E[X] = np, Var[X] = np(1-p).

clearvars; close all
rng(0)

n = 20;         % number of Bernoulli draws in one experiment
p = 0.3;        % success probability
N = 1e5;        % number of experiments

b = rand(n, N) < p;     % Bernoulli draws, one experiment per column
x = sum(b, 1);          % number of successes in each experiment

k = 0:n;
p_emp = sum(x == k', 2)'/N;
p_theory = arrayfun(@(kk) nchoosek(n, kk), k) .* p.^k .* (1-p).^(n-k);

figure
bar(k, p_emp)
hold on
stem(k, p_theory, 'r', 'LineWidth', 1.5)
hold off
xlabel('k')
ylabel('p_X[k]')
legend('simulation', 'theory, Bin(n,p)')
grid on

fprintf('E[X]:   theory = %.3f, simulation = %.3f\n', n*p, mean(x))
fprintf('Var[X]: theory = %.3f, simulation = %.3f\n', n*p*(1-p), var(x))