function q=qskt(p,location,scale,shape,df,tol)
%dskt pskt qskt rskt
%Skew-t Distribution
% 
%DESCRIPTION
% 
%Density function, distribution function, quantiles and random number
%generation for the skew-t (S-t) distribution.
% 
%USAGE
% 
%dskt(x, location, scale, shape, df)
%pskt(q, location, scale, shape, df)
%qskt(p, location, scale, shape,  df, tol)
%rskt(n, location, scale, shape, df)
% 
%REQUIRED ARGUMENTS
% 
%x	vector of quantiles. Missing values (NaN) are allowed.
%q	vector of quantiles. Missing values (NaN) are allowed.
%p	vector of probabilities. Missing values (NaN) are allowed.
%n  scalar value for the sample size  
%
%OPTIONAL ARGUMENTS
% 
%location vector of location parameters (default is 0).
%scale	  vector of (positive) scale parameters (default is 1).
%shape	  vector of shape parameters. With 'pskt' and 'qskt', it must 
%	      be of length 1 (default is 0).
%n	      sample size (default is 1).
%df       scalar value for the degrees of freedom (default is positive 
%         infinity (Inf) which corresponds to the skew normal distribution)
%tol	  a scalar value which regulates the accuracy of the result.
%         (default is 1e-8) 

%VALUE
% 
%density (dskt), probability (pskt), quantile (qskt) or  random sample (rskt)
%from the skew-t distribution with given location, scale, shape and degrees
%of freedom parameters.
% 
%BACKGROUND
% 
%The family of skew-t distributions is an extension of the Student's t family, 
%via the introduction of a shape parameter which regulates skewness; when shape=0, 
%the skew-t distribution reduces to the usual Student's t distribution. When df=Inf,
%it reduces to the skew normal distribution. A multivariate version of the distribution 
%exists. See the reference below for additional information.
% 
%DETAILS
% 
%pskt make use of function T_Owen
%qskt make use of functions: (1) dskt (2)pskt (3)skt_cumulants 
% 
%REFERENCES
%
%Azzalini, A. and Capitanio,A. (2003). Distributions generated by perturbation of 
%symmetry with emphasis on a multivariate skew-t distribution. J.Roy.Statist.Soc.
%B 65, 367-389.
%
%Azzalini, A. (1985). A class of distributions which includes the normal
%ones. Scand. J. Statist. 12, 171-178.
% 
%Azzalini, A. and Dalla Valle, A. (1996). The multivariate skew-normal
%distribution. Biometrika 83, 715-726.
% 
%SEE ALSO
% 
%T_Owen, 
% 
%EXAMPLES
% 
%pdf = dskt(-4:0.1:4,[],[],3,5)
%rnd = rskt(100, 5, 2, -5, 8) 
%q   = qskt([0.25 0.5 0.75],0,1,3,5,[])
%pst = pskt(q, [], [], 3, 5)        #should give back pst=[0.25 0.5 0.75]
%p   = [nan 2 -1 0 0.25 1 nan 0.75 1 0.5]
%q   = qskt(p,0,1,2,5,[])

if nargin<6
    tol=1e-08;
end;
if nargin<5
    df=Inf;
end;
if nargin<4
    shape=0;
end;
if nargin<3
    scale=1;
end;
if nargin<2
    location=0;
end;
if nargin<1
    error('argument "x" is missing, with no default');
end;

%Set the default parameter to let the program work with Nan
%Nan=(Not-available-number)
if isnan(tol)
    tol=1e-08;
end;
if isnan(df)
    df=Inf;
end;
if isnan(shape)
    shape=0;
end;
if isnan(scale)
    scale=1;
end;
if isnan(location)
    location=0;
end;

%Set the default parameter values for the empty values: []
if isempty(tol)
    tol=1e-08;
end;
if isempty(df)
    df=Inf;
end;
if isempty(shape)
    shape=0;
end;
if isempty(scale)
    scale=1;
end;
if isempty(location)
    location=0;
end;
if isempty(p)
    error('You need to provide a vector of probabilities');
end;

%Define constraints on the scale parameter
if scale<=0
    error('The scale parameter must have a positive nonzero value');
end;

if df==Inf
    q=qsn(p,location,scale,shape);
    elseif df~=Inf;
    %Define constraints on the DF parameter
        if df<1 || df~=round(df)
        error('degrees of freedom must be an integer positive number or infinity');
        end;

        %begin program computations
        max_q=sqrt(finv(p,1,df));
        min_q=-sqrt(finv(1-p,1,df));

        if shape==Inf
            q=location+scale*max_q;
        elseif shape==-Inf
            q=location+scale*min_q;
        elseif shape~=Inf || shape~=-Inf
            %Set the values of the vector of probabilities such that it admits nan's
            a=isnan(p);
            ind_nan=find(a==1);
            ind_p1=find(p==1);
            ind_p0=find(p==0);

            ind_p_L_0=find(p<0);
            ind_p_H_1=find(p>1);

            p(ind_nan)=0.5;
            p(ind_p1)=0.5;
            p(ind_p0)=0.5;
            p(ind_p_L_0)=0.5;
            p(ind_p_H_1)=0.5;

            cum=skt_cumulants(0,1,shape,max(df,5),4);
            g1=cum(3)/cum(2)^(3/2);
            g2=cum(4)/cum(2)^2;
            x=norminv(p);
            x=(x+(x.^2-1).*g1./6+x.*(x.^2-3).*g2./24-x.*(2.*x.^2-5).*g1.^2./36);
            x=cum(1)+sqrt(cum(2))*x;
            max_err=1;
            while (max_err>tol)
                x1=x-(pskt(x,0,1,shape,df)-p)./dskt(x,0,1,shape,df);
                x1=min(x1,max_q);
                x1=max(x1,min_q);
                max_err=max(abs(x1-x)./(1+abs(x)));
                x=x1;
            end;

            x(ind_nan)=nan;
            x(ind_p_L_0)=nan;
            x(ind_p_H_1)=nan;
            x(ind_p1)=Inf;
            x(ind_p0)=-Inf;

            q=location+scale*x;
        end;
end;


