function [r,type,coefs] = johnsrnd(quantiles,varargin)
%JOHNSRND Random arrays from the Johnson system of distributions.
%   R = JOHNSRND(QUANTILES,M,N) returns an M-by-N matrix of random numbers
%   drawn from the distribution in the Johnson system that satisfies the
%   quantile specification given by QUANTILES.  QUANTILES is a 4-element
%   vector of quantiles for the desired distribution that correspond to the
%   standard normal quantiles [-1.5 -0.5 0.5 1.5].  In other words, you
%   specify a distribution from which to draw random values by designating 
%   quantiles that correspond to the cumulative probabilities [0.067 0.309
%   0.691 0.933]. QUANTILES may also be a 2x4 matrix whose first row contains
%   4 standard normal quantiles, and whose second row contains the
%   corresponding quantiles of the desired distribution.  The standard normal
%   quantiles must be spaced evenly.
%
%   Note: Because R is a random sample, its sample quantiles will typically
%   differ somewhat from the specified distribution quantiles.
%
%   R = JOHNSRND(QUANTILES) returns a scalar value.
%   R = JOHNSRND(QUANTILES,M,N,...) or R = JOHNSRND(QUANTILES,[M,N,...])
%   returns an M-by-N-by-... array.
%
%   [R,TYPE] = JOHNSRND(...) returns the type of the specified distribution
%   within the Johnson system.  TYPE is 'SN', 'SL', 'SB', or 'SU'.  Set M
%   and N to zero to identify the distribution type without generating any
%   random values.
%
%   The four distribution types in the Johnson system correspond to the
%   following transformations of a normal random variate:
%
%      Type SN: Identity transformation (normal distribution)
%      Type SL: Exponential transformation (lognormal distribution)
%      Type SB: Logistic transformation (bounded)
%      Type SU: Hyperbolic sine transformation (unbounded)
%
%   [R,TYPE,COEFS] = JOHNSRND(...) returns coefficients of the transformation
%   that defines the distribution.  COEFS is [GAMMA, ETA, EPSILON, LAMBDA].  If
%   Z is a standard normal random variable and H is one of the transformations
%   defined above, then R = LAMBDA * H((Z - GAMMA) / ETA) + EPSILON is a random
%   variate from the distribution type corresponding to H.
%
%   Examples
%      % Generate random values that with longer tails than a standard normal
%      r = johnsrnd([-1.7 -.5 .5 1.7],1000,1);
%      qqplot(r);
%
%      % Generate random values skewed to the right
%      r = johnsrnd([-1.3 -.5 .5 1.7],1000,1);
%      qqplot(r)
%
%      % Generate random values that match some sample data well in the
%      % right-hand tail
%      load carbig;
%      qnorm = [.5 1 1.5 2];
%      q = quantile(Acceleration, normcdf(qnorm));
%      r = johnsrnd([qnorm; q],1000,1);
%      [q; quantile(r, normcdf(qnorm))]
%
%      % Determine the distribution type and the coefficients
%      [r,type,coefs] = johnsrnd([qnorm; q],0);  % returns [] for r
%
%   See also RANDOM, PEARSRND.

%   JOHNSRND uses transformations of standard normal random variates.

%   References:
%      [1] Johnson, N.L., S. Kotz, and N. Balakrishnan (1994) Continuous
%          Univariate Distributions, Volume 1,  Wiley-Interscience.
%      [2] Johnson, N.L. (1949) Systems of Frequency Curves Generated by
%          Methods of Translation, Biometrika 36:149-176.
%      [3] Slifker, J.F. and S.S. Shapiro (1980) "The Johnson System:
%          Selection and Parameter Estimation", Technometrics 22(2):239-246.
%      [4] Wheeler, R.E. (1980) "Quantile estimators of Johnson curve
%          parameters", Biometrika 67(3):725�728.

%   Copyright 2005 The MathWorks, Inc. 
%   $Revision: 1.1.6.1 $  $Date: 2005/11/18 14:28:08 $

if nargin < 1
    error('stats:johnsrnd:TooFewInputs','Requires at least one input argument.');
elseif isvector(quantiles)
    if numel(quantiles) ~= 4
        error('stats:johnsrnd:BadQuantiles','QUANTILES must be a 4-element vector or a 2x4 matrix.');
    elseif ~all(diff(quantiles) > 0)
        error('stats:johnsrnd:BadQuantiles','QUANTILES must contain montonically increaing values.');
    end
    z0 = 0;
    dz = 0.5;
elseif isequal(size(quantiles),[2,4])
    z = quantiles(1,:);
    diffz = diff(z);
    quantiles = quantiles(2,:);
    if ~all(diff(quantiles) > 0)
        error('stats:johnsrnd:BadQuantiles','The second row of QUANTILES must contain montonically increasing values.');
    elseif ~all(diffz>0)
        error('stats:johnsrnd:BadQuantiles','The first row of QUANTILES must contain monotonically increasing values.');
    elseif ~(range(diffz) < eps(max(diffz)))
        error('stats:johnsrnd:BadQuantiles','The first row of QUANTILES must contain evenly-spaced normal quantiles.');
    end
    z0 = mean(z);
    dz = .5*mean(diffz);
else
    error('stats:johnsrnd:BadQuantiles','QUANTILES must be a 4-element vector or a 2x4 matrix.');
end

[err, sizeOut] = statsizechk(1,1,varargin{:});
if err > 0
    error('stats:johnsrnd:InputSizeMismatch','Size information is inconsistent.');
end

% Use Slifker and Shapiro's algorithm to determine which of the four cases
m = quantiles(4) - quantiles(3);
n = quantiles(2) - quantiles(1);
p = quantiles(3) - quantiles(2);
q0 = .5*(quantiles(2)+quantiles(3));

mp = m./p;
np = n./p;

tol = 1e-4;

% Compute the two sets of location-scale parameters according to which case
% we're in. Then generate standard normals and transform.
%
% The bounded case: Z = gamma + eta*log((R-epsilon)/(lambda-y+epsilon))
if mp.*np < 1 - tol
    type = 'SB'; % bounded
    
    pm = p./m;
    pn = p./n;
    eta = dz ./ acosh(.5*sqrt((1+pm).*(1+pn)));
    gamma = z0 + eta .* asinh((pn-pm).*sqrt((1+pm).*(1+pn)-4) ./ (2*(pm.*pn-1)));
    lambda = p .* sqrt(((1+pm).*(1+pn) - 2).^2 - 4) ./ (pm.*pn-1);
    epsilon = q0 - .5*lambda + p.*(pn-pm) ./ (2*(pm.*pn-1));
    
    w = (randn(sizeOut) - gamma) / eta;
    v = exp(-abs(w));
    signw = sign(w) + (w == 0);
    r = lambda .* (signw.*abs((1-v)./(1+v)) + 1)./2 + epsilon;
    
% The unbounded case: Z = gamma + eta*asinh((R-epsilon)/lambda)
elseif mp.*np > 1 + tol
    type = 'SU'; % unbounded

    eta = 2*dz ./ acosh(.5*(mp+np));
    gamma = z0 + eta .* asinh((np - mp) ./ (2*sqrt(mp.*np-1)));
    lambda = 2*p .* sqrt(mp.*np-1) ./ ((mp+np-2).*sqrt(mp+np+2));
    epsilon = q0 + p .* (np-mp)./(2*(mp+np-2));

    w = exp((randn(sizeOut) - gamma) ./ eta);
    r = .5*lambda.*(w - 1./w) + epsilon;
    
else % abs(mp.*np - 1) < tol
    % Lognormal case: Z = gamma + eta*log(R-epsilon)
    if abs(mp - 1) > tol
        type = 'SL'; % lognormal

        % When mp < 1, gamma/eta will have -pi imaginary part.  After
        % exponentiating, that's equivalent to lambda == -1 instead of 1.
        eta = 2*dz ./ log(mp);
        gamma = z0 + eta .* log(abs(mp-1) ./ (p.*sqrt(mp)));
        lambda = sign(mp-1);
        epsilon = q0 - .5*p .* (mp+1)./(mp-1);

        r = lambda.*exp((randn(sizeOut) - gamma) ./ eta) + epsilon;
        
    % Limiting normal case: Z = (R-epsilon)/lambda
    else
        type = 'SN'; % normal

        eta =  2*dz./m;
        gamma = z0 - q0.*eta;
        lambda = 1;
        epsilon = 0;
        
        r = (randn(sizeOut) - gamma) ./ eta;
    end
end

coefs = [gamma eta epsilon lambda];