function x = generate_nums_ber(n, p ,m)
% GENERATE_NUM_BER - generate m outcomes using n bernulli trials

%Create vector of size m, go through loop m times
x=zeros(1,m); % 1 x m vector of zeros for storing outcomes
for j=1:m
  successful_trials = 0;
  for i = 1:n
      u = rand; %between [0,1)
      if(u <= p)
          successful_trials= successful_trials+1;
      end
  end
  x(j)= successful_trials;
end

% alternative to loop above: x=(rand(1,m)>=1-p);

% plot normalised histogram of outcomes versus theoretical pmf

close all;  % close all previous plots

alphabet_x=zeros(1,n);        % 1 x n vector with alphabet of r.v.
for i = 1:n+1
    alphabet_x(i) = i-1;
end

hx=histc(x,alphabet_x);  % compute histogram of outcomes x on alphabet
stem(alphabet_x,hx/m);   % plot normalised histogram on alphabet
hold on;                 % keep this figure when plotting next
px=zeros(1, n+1);              % theoretical pmf of Ber(p)
for i = 0:n
    px(i+1) = binopdf(i, n, p);
end
stem(alphabet_x,px,'rs'); % plot pmf

axis([-1 n+1 0 1]);                  % widen axis of figure
legend('empirical','theoretical'); % add legend to figure