function r=rskt(n,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.
% 
%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 'psn' and 'qsn', it must 
%	      be of length 1 (default is 0).
%n	      sample size (default is 1).
%tol	  a scalar value which regulates the accuracy of the result.
%         (default is 1e-8) 
%df       scalar value for the degrees of freedom (default is positive 
%         infinity which corresponds to the skew normal distribution)
%
%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
%rskt make use of function rsn
%
%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 <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 n is missing');
end;


if isnan(df);
    df=Inf;
end;

if isnan(location);
    location=0;
end;

if isnan(scale);
    scale=1;
end;

if isnan(shape);
    shape=0;
end;

if isnan(n);
    n=1;
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(n)
    n=1;
end;

%Define constraints on the scale parameter
if scale<=0
    error('The scale parameter must have a positive nonzero value');
end;
%Define constraints on the sample size
if n<1
    error('The default sample size must be equal to 1');
end;

if round(n)~=n
    error('The default sample size must be an integer number');
end;

%Display the parameters that have been used to generate the random numbers
A = {'Parameters',' ';'Location',location;'Scale',scale;'Shape',shape;'Degrees of Freedom',df};
%display(A);

if df==Inf;
    %Creates a (1*n) matrix with random numbers from a skew normal distribution
    %where the location parameter is set equal to zero
    z=rsn(n,0,scale,shape);
    r=z+location;
elseif df~=Inf;
    %Creates a (1*n) matrix with random numbers from a skew normal distribution
    %where the location parameter is set equal to zero
    z=rsn(n,0,scale,shape);
    v=chi2rnd(df,1,n)./df;
    y=z./sqrt(v)+location;
    r=y;
end
end


function y = rsn(n,xi,omega,alpha)
    delta = alpha/sqrt(1+alpha^2);
    tn    = randn(2,n);
    chi   = abs(tn(1,:));
    nrv   = tn(2,:);
    z     = delta*chi+sqrt(1-delta^2)*nrv;
    y     = xi+omega*z;
end



% rst <- function (n=1, xi = 0, omega = 1, alpha = 0, nu=Inf, dp=NULL)
% { 
%     if(!is.null(dp)) {
%      if(!missing(alpha)) 
%         stop("You cannot set both 'dp' and the component parameters")
%      xi <- dp[1]
%      omega <- dp[2]
%      alpha <- dp[3]
%      nu <- dp[4]
%     }
%   # rs <<- .Random.seed
%   z <- rsn(n, 0, omega, alpha)
%   if(nu < Inf) {  
%     # .Random.seed <<- rs  
%     v <- rchisq(n,nu)/nu
%     y <- z/sqrt(v) + xi
%     }
%     else y <- z + xi
%   attr(y, "family") <- "ST"
%   attr(y, "parameters") <- c(xi, omega, alpha, nu)
%   return(y)
% }



% rsn <- function (n = 1, xi = 0, omega = 1, alpha = 0, tau = 0, dp = NULL) 
% {# since version 1.6-2 (2020): use transformation/additive method throughout
%   if (!is.null(dp)) {
%     if (!missing(alpha)) 
%         stop("You cannot set both 'dp' and the component parameters")
%     xi <- dp[1]
%     omega <- dp[2]
%     alpha <- dp[3]
%     tau <- if (length(dp) > 3) dp[4] else 0
%   }  
%   delta <- alpha/sqrt(1 + alpha^2)
%   if(tau == 0) {
%     tn <- matrix(rnorm(2*n), 2, n,  byrow = FALSE)
%     chi <- c(abs(tn[1,]))
%     nrv <- c(tn[2,])
%     z <- delta * chi + sqrt(1 - delta^2) * nrv
%   } else {
%     # rs <<- .Random.seed
%     truncN <- qnorm(runif(n, min= pnorm(-tau), max=1))
%     # .Random.seed <<- rs
%     z <- delta * truncN + sqrt(1-delta^2) * rnorm(n)
%   }
%   y <- as.vector(xi + omega * z)
%   attr(y, "family") <- "SN"
%   attr(y, "parameters") <- c(xi, omega, alpha, tau)
%   return(y)
% }

