function x = quadratic_rnd(a,b,r,c)
% x = quadratic_rnd( A, B, C, r,c)
%
% generate a matrix of a random variable distribued according to the 
% quadratic distribution with parameters  a and b
%
%  INPUT         DESCRIPTION
%   a            lower bound of X                             1 x 1   or  r x c
%   b            upper bound of X                             1 x 1   or  r x c
%   r            number of rows in x                                      1 x 1
%                sample of standard normal variables, possibly correlated r x c
%   c            number of cols in x                                      1 x 1
%
%  OUTPUT
%   x            sample of the extreme type I distribution                r x c
%
%  EXAMPLE 
%    To plot the PDF of a single quadratically distributed random variable, ... 
%    n=10000; x=quadratic_rnd(2,3,1,n); hist(x,n/100,1); avg_x = sum(x)/n

 u = rand(r,c);

%                           3      2                    2        2      3
%     f := - U*(a-b)^3/6 + a  - 3 a  b + 6 a b x - 3 a x  - 3 b x  + 2 x  

for ii = 1:r
  for jj = 1:c
    % coefficients of cubic equation
    C = [ 2 , -3*(a+b) , 6*a*b , a^3 - 3*a^2*b - u(ii,jj)*(a-b)^3 ];
    % roots of cubic equation
    R = roots(C);
    % root of cubic equation between a and b
    x(ii,jj) = R(find( (a<R) & (R<b)));
  end
end

% ----------------------------------------------------------  quadratic_rnd
% HP Gavin 2023-02-14
