function x = generate_geometric_bernoulli(p, m)
  % generate geometric RV X with parameter p, m times using Bernoulli Trials

  emperical = zeros(1,m);                       % allocate array for emperical results of m simulations
  for i = 1:m                                   % iterate over m simulations
    failures = 0;                               % count number of failures before a successful trial is achieved
    u = rand;                                   % generate random nuumber from 0-1
    while (u > p)                               % Loop until succesful trial
      failures++;                               % count each unsuccessful trial
      u = rand;                                 % generate new random number for next trial
    endwhile
    emperical(i) = failures + 1;                % add the number of trials taken for a successful trial
  end

  close all;                                    % close any existing graphs

  max_x = max(emperical);
  x_values = [1 : max_x];                       % array of x-axis values
  hist(emperical, x_values, "facecolor", "r");  % plot empirical data
  xlim([-0.5 (max_x + 0.5)]);                   % set x-axis to allow for histogram bar widths

end