function rnd = discrete_rnd (v, p, varargin)
% discrete_rnd (v, p)
% discrete_rnd (v, p, r)
% discrete_rnd (v, p, r, c, ...)
% discrete_rnd (v, p, [sz])
% Return a matrix of random samples from the univariate distribution which
% assumes the values in v with probabilities p.
%
% When called with a single size argument, return a square matrix with
% the dimension specified.  When called with more than one scalar argument the
% first two arguments are taken as the number of rows and columns and any
% further arguments specify additional matrix dimensions.  The size may also
% be specified with a vector of dimensions sz.
%
% If no size arguments are given then the result matrix is the common size of
% v and p.

% Author: KH <Kurt.Hornik@wu-wien.ac.at>
% Description: Random deviates from a discrete distribution
% modified by Andreas Stahel to work on Matlab too

  if (nargin < 2)
    print_usage ();
  end%if

  if ( ~isvector (v))
    error ('discrete_rnd: V must be a vector');
  elseif ( ~isvector (p) || (length (p) ~= length (v)))
    error ('discrete_rnd: P must be a vector with length (V) elements');
  elseif (any (isnan (p)))
    error ('discrete_rnd: P must not have any NaN elements');
  elseif ( ~(all (p >= 0) && any (p)))
    error ('discrete_rnd: P must be a nonzero, non-negative vector');
  end%if

  if (nargin == 2)
    sz = size (v);
  elseif (nargin == 3)
    if (isscalar (varargin{1}) && varargin{1} >= 0)
      sz = [varargin{1}, varargin{1}];
    elseif (isrow (varargin{1}) && all (varargin{1} >= 0))
      sz = varargin{1};
    else
      error ('discrete_rnd: dimension vector must be row vector of non-negative integers');
    end%if
  elseif (nargin > 3)
    if (any (cellfun (@(x) (~isscalar (x) || x < 0), varargin)))
      error ('discrete_rnd: dimensions must be non-negative integers');
    end%if
    sz = [varargin{:}];
  end%if

  if (~isvector (v))
    error ('discrete_rnd: V must be a vector');
  end%if
  if (~isvector (p) || (length (p) ~= length (v)))
    error ('discrete_rnd: P must be a vector with length (V) elements');
  end%if
  if (any (isnan (p)))
    error ('discrete_rnd: P must not have any NaN elements');
  end%if
  if (~(all (p >= 0) && any (p)))
    error ('discrete_rnd: P must be a nonzero, non-negative vector');
  end%if
  
  if (nargin == 2)
    sz = size (v);
  end%if
  if (nargin == 3)
    if (isscalar (varargin{1}) && varargin{1} >= 0)
      sz = [varargin{1}, varargin{1}]
    elseif (isrow (varargin{1}) && all (varargin{1} >= 0))
      sz = varargin{1};
    else
      error ('discrete_rnd: dimension vector must be row vector of non-negative integers');
    end%if (isrow
  end%if (nargin
  if (nargin > 3)
    if (any (cellfun (@(x) (~isscalar (x) || x < 0), varargin)))
      error ('discrete_rnd: dimensions must be non-negative integers');
    else
      sz = [varargin{:}];
    end%if (any
  end%if (nargin
  
  rnd = rand(sz); rnd = rnd(:);
  %rnd = v(sum(rnd(:) >= cumsum(p(1:end-1))/sum(p),2)+1);
  res = [];
  for lim = cumsum(p(1:end-1)/sum(p))
    res = [res rnd > lim];
  end%for

  rnd = v(sum(res,2)+1);
  rnd = reshape (rnd, sz);
end%function
