function mbetarv = genbeta(nrows,ncols,palpha,pbeta) ;
% GENBETA, Generates Beta(palpha,pbeta) psuedo random variables
% Inputs:
%      nrows   - number of rows in output matrix
%      ncols   - number of columns in output matrix
%      palpha  - scalar parameter alpha, assumed >= 1
%      pbeta   - scalar parameter beta, assumed >= 1
% Output:
%     mbetarv  - matrix of beta pseudo random realizations
% Note:
%     Control the seed with rand('seed', ... )

ntot = nrows * ncols ;
          %  TOTAL NUMBER OF OBSERVATIONS REQUIRED

%  ---------------- GET MAXIMUM VALUE OF DENSITY ----------------
if palpha == 1 ;
          %  THEN MAX IS AT LEFT ENDPOINT
  fmax = pbeta ;
          %  SINCE MAX VALUE OF DENSITY IS
          %                      gamma(pbeta + 1) / gamma(pbeta)
elseif pbeta == 1 ;
          %  THEN MAX IS AT RIGHT ENDPOINT
  fmax = palpha ;
          %  SINCE MAX VALUE OF DENSITY IS
          %                      gamma(palpha + 1) / gamma(alpha)
else ;
          %  THEN MAX IS IN INTERIOR OF (0,1)
  xmax = (palpha - 1) / (palpha + pbeta - 2) ;
          %  LOCATION OF MAX
  fmax = fbeta(xmax,palpha,pbeta) ;
          %  MAX VALUE OF BETA DENSITY
end ;

% ----------------- GENERATE FIRST BLOCK OF DATA --------------------
%                TO GET ESTIMATE OF REJECTION PROB

rvpb = max([ntot; 100]) ;
          %  RV'S PER BLOCK
          %  AT LEAST 100 FOR DECENT EST. OF REJ. PROB.
rvpb = min([rvpb; 1000]) ;
          %  LESS THAN 1000 SHOULD BE ENOUGH TO EST PROB.
udata = rand(rvpb,2) ;
          %  UNIF'S ON (0,1)x(0,1)
udata(:,2) = udata(:,2) * fmax ;
          %  UNIF'S ON (0,1)x(0,MAX(F))
flag = (udata(:,2) < fbeta(udata(:,1),palpha,pbeta)) ;
          %  FLAG (WITH 1, 0 OTHERWISE) ROWS WHERE WILL NOT REJECT
eprob = mean(flag) ;
          %  ESTIMATED CHANCE OF NONREJECTION
rb = udata(flag,1) ;
          %  START OF VECTOR OF VECTOR OF RANDOM BETA'S
          %  RECALL NONREJECTED X'S HAVE BETA DISTRIBUTION
nobs = length(rb) ;

if eprob == 0 ;
  eprob = 1 / (rvpb + 1) ;
          %  BEST GUESS ABOUT REJECTION PROB WHEN DON'T SEE ANY
end ;

disp(['Finished 1st block, have ' num2str(nobs) ' obs''s so far']) ;


% ------------------- LOOP THROUGH NEXT BLOCKS -------------------
while nobs < ntot ;
  rvpb = roots([eprob; -2*sqrt(eprob*(1-eprob)); nobs - ntot]) ;
          %  WANT rvpb AT EACH POINT TO BE 95% SURE OF PROVIDING
          %  ENOUGH OBSERVATIONS, I.E. NEED TO SOLVE FOR n, THE
          %  QUADARTIC (IN sqrt(n)):   
          %          np - 2sqrt(np(1-p)) = ntot - nobs
  rvpb = (max(real(rvpb)))^2 ;
          %  SQUARE THE BIGGER ROOT OF THE QUADRATIC
          %  USE real SINCE CAN HAVE SLIGHT COMPLEX PART DUE TO 
          %  ROUND OFF ERROR
  rvpb = ceil(rvpb) ;
          %  ROUND UP TO THE NEXT INTEGER
  rvpb = max([10; min([rvpb; 10000])]) ;
          %  PERHAPS INEFFICIENT TO DO LESS THAN 10 OR MORE THAN 10,000

  %  NOW USE LINES FROM ABOVE TO ADD TO VECTOR OF BETA'S
  udata = rand(rvpb,2) ;
  udata(:,2) = udata(:,2) * fmax ;
  flag = (udata(:,2) < fbeta(udata(:,1),palpha,pbeta)) ;
  rb = [rb; udata(flag,1)] ;
  nobs = length(rb) ;

  disp(['Finished next block, have ' num2str(nobs) ' obs''s so far']) ;
end ;

rb = rb(1:ntot) ;
          %  ONLY NEED FIRST ntot OF THEM
mbetarv = reshape(rb,nrows,ncols) ;
          %  PUT INTO FINAL MATRIX FORM
