function x = corr_beta_rnd( a, b, q, p, R, N )
% x = corr_beta_rnd( a, b, q, p, R, N )
%
% generate a sample of m observations of n beta random variables to the 
% with parameters  a, b, p and q
% the standard normal deviates of X are correlated with matrix R
%
%  INPUT         DESCRIPTION
%   a            lower bound of x                                         n x 1
%   b            upper bound of x                                         n x 1
%   q            exponent                                                 n x 1
%   p            exponent                                                 n x 1
%   R            correlation matrix of standard normal deviates           n x n
%   N            number of observations                                   1 x 1
%
%    Values of q and p must be greater than zero.  
%
%  OUTPUT
%   x            sample of m observations of n correlated betaormal r.v's n x m
%
%    If X is a betaormal random variable, then the log of X is normally 
%    distributed with a mean of  log(medX) and a variance of log(1+covX^2)
%
%  EXAMPLE 
%    To plot the PDF of a single betaormal random variable, ... 
%    n=10000; x=beta_rnd(1,0.5,1,n); hist(x,n/100,1); avg_x = sum(x)/n
 
%    Reference: http://en.wikipedia.org/wiki/Beta_distribution

 if nargin ~= 6
    help corr_logn_rnd
    return
 end

 n = length(a);      % number of random variables

 if any(diag(R)~=1), error(' corr_logn_rnd: diagonal of R must equal 1'); end
 if any(R>1)  error(' corr_logn_rnd: R values cannot be greater than 1'); end
 if any(R<-1) error(' corr_logn_rnd: R values cannot be less than 1');   end
 if size(b,1)~=n error(sprintf(' b must have %d rows, zero columns',n)); end
 if size(q,1)~=n error(sprintf(' q must have %d rows, zero columns',n)); end
 if size(p,1)~=n error(sprintf(' q must have %d rows, zero columns',n)); end
 if size(R,1)~=n error(sprintf(' R must have %d rows, %d columns',n,n)); end
 if size(R,2)~=n error(sprintf(' R must have %d rows, %d columns',n,n)); end
 if any(q<0) error(' q must be positve'); end
 if any(p<0) error(' p must be positve'); end

% eigenvalues and eigenvectors of the correlation matrix R
 [eVec,eVal] = eig(R);
 if any(eVal<0), error(' corr_beta_rnd:  R must be positive definite'); end

 % generate N observations of n uncorrelated beta variables

 x = NaN(n,N);
 Y = norm_cdf(eVec* sqrt(eVal) * randn(n,N));
 for ii = 1:n
   x(ii,:) = a(ii) + (b(ii)-a(ii))*betaincinv( Y(ii,:) , q(ii) , p(ii) );
 end

% ------------------------------------------------------------- corr_beta_rnd
% HP Gavin  2024-03-07
