function P=pmskt(x, xi, Omega, alpha, df)
%
%Description
%Probability density function, distribution function and random number 
%generation for the multivariate skew-t (MST) distribution. 
%
%Usage
% dmst(x, xi=rep(0,length(alpha)), Omega, alpha, df=Inf, log=FALSE)
% dmst(x, dp=, log=FALSE)
% pmskt(x, xi=rep(0,length(alpha)), Omega, alpha, df=Inf, ...)
% pmskt(x, dp=, ...)
% rmst(n=1, xi=rep(0,length(alpha)), Omega, alpha, df=Inf)
% rmst(n=1, dp=)
%
%Arguments
% x      for dmsn, this is either a vector of length d, where d=length(alpha), 
%        or a matrix with d columns, giving the coordinates of the point(s) where 
%        the density must be evaluated; for pmsn, only a vector of length d is allowed.  
% xi     a numeric vector of lenght d, or a matrix with d columns, representing 
%        the location parameter of the distribution. If xi is a matrix, its dimensions 
%        must agree with those of x.  
% Omega  a positive-definite covariance matrix of dimension (d,d).  
% alpha  a numeric vector which regulates the shape of the density.  
% df     degrees of freedom (scalar); default is df=Inf which corresponds to the 
%        multivariate skew-normal distribution 
% n      a numeric value which represents the number of random vectors to be drawn.  
%
%Details
% The positive-definiteness of Omega is not tested for efficiency reasons. 
% Function pmst requires pmt from package mnormt; the accuracy of its computation can 
% be controlled via use of ... 
%
%Value
% A vector of density values (dmst) or a single probability (pmst) or a matrix of random 
% points (rmst).
%
%Background
% The family of multivariate skew-t distributions is an extension of the multivariate 
% 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 t distribution. When df=Inf 
% the distribution reduces to the multivariate skew-normal one; see dmsn. See the 
% reference below for additional information. 
%
%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. 
%
%See Also
%dst, dmsn, dmt 

if nargin<5; 
    df=Inf;
end;

if nargin<2; 
    xi=repmat(0,1,length(x));
end;

if nargin<1;
    error('Argument x is missing');
end;

if isnan(xi); 
    xi=repmat(0,1,length(x));
end;

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


d=length(alpha);
omega=sqrt(diag(Omega))';

%Get the correlation matrix from the covariance matrix Omega
varm=inv(sqrt(diag(diag(Omega))));
Corr_M=varm*Omega*varm;

Ocor=Corr_M;
O_alpha=(Ocor*alpha')';
delta=O_alpha/sqrt(1+sum(alpha.*O_alpha));
%equivalent of rbind
q=[1 -delta];
%equivalent of cbind
w=[-delta' Ocor];
%equivalent of rbind(...)
Obig=[q' w'];


y=size(x);
n=y(1,1);

Z=repmat(xi',n,1);
Z1=reshape(Z,n,d);
%Equivalent of x-xi
Z2=x-Z1;
Z3=repmat(omega',n,1);
Z4=reshape(Z3,n,d);
%Equivalent of (x-xi)/omega
Z5=Z2./Z4;
Z6=reshape(Z5,1,n*d);
x=[0 Z6];

if df==Inf
    mu=zeros(1,d+1);
    P=2*mvncdf(x,mu,Obig);
else
    P=2*mvtcdf(x,Obig,df);
end

    



%Examples
% x <- seq(-4,4,length=15)
% xi <- c(0.5, -1)
% Omega <- diag(2)
% Omega[2,1] <- Omega[1,2] <- 0.5
% alpha <- c(2,2)
% pdf <- dmst(cbind(x,2*x-1), xi, Omega, alpha, df=5)
% rnd <- rmst(10,  xi, Omega, alpha, 6)
% p1 <- pmst(c(2,1), xi, Omega, alpha, df=5)
% p2 <- pmst(c(2,1), xi, Omega, alpha, df=5, abseps=1e-12, maxpts=10000)





