function d=dskt(x,location,scale,shape,df)
%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 which corresponds to the skew normal distribution)
%tol	  a scal 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
% 
%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
%
%a   =dskt(1)                     #should give back 0.2420
%b   =dskt(1,0,1,0,Inf)           #should give back 0.2420
%c   =dskt(1,[],[],[],[])         #should give back 0.2420
%d   =dskt()                      #error: Argument x is missing
%e   =dskt([])                    #error: x cannot be an empty vector
%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<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');
end;

if isnan(location); 
    location=0;
end;
if isnan(scale);
    scale=1;
end;
if isnan(shape);
    shape=0;
end;
if isnan(df);
    df=Inf;
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(x)
    error('x cannot be an empty vector');
end;

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

%Define the standardized random variable
if df==Inf;
    z=(x-location)./scale;
    d=2.*normpdf(z).*normcdf( shape.*z)./scale;
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;
    z=(x-location)./scale;
    d=2.*tpdf(z,df).*tcdf(shape.*z.*sqrt((1+df)./(z.^2+df)),df+1)./scale;
end;

