[R] Multivariate skew-t cdf

[Spencer Graves]

	  You want to evaluate skewed t probabilities in how many dimensions? 
If 27 is your maximum, the problem won't be as difficult as if 27 is 
your minimum.  Also, do you want to compute the multivariate cumulative 
probability function for arbitrary points, location, covariance, shape 
and degrees of freedom?  Or are you really only interested in certain 
special case(s)?  If you have a simpler special case, it might be easier 
to get a solution.

	  I was able to replicate your result, even when I reduced the 
dimensionality down to 20;  with 19 dimensions, the function seemed to 
return a sensible answer.  If it were my problem, I might first make a 
local copy of the pmst function and modify it to use the mvtnorm package 
in place of mnormt.  That might get you answers with somewhat higher 
dimensionality, though it might not be adequate -- and I wouldn't trust 
the numbers I got without serious independent testing.  I'd try to think 
how I could modify the skewness so I could check the accuracy that way.

	  Have you studied the reference mentioned in the "dmst" help page, and 
reviewed some of their sources?  Computing multivariate probabilities 
like this is still a research project, I believe.  In this regard, I 
found the following two books helpful:

	  * Evans and Schwarz (2000) Approximating Integrals via Monte Carlo 
and Deterministic Methods (Oxford)

	  * Kythe and Schaeferkotter (2005) Handbook of Computational Methods 
for Integration (Chapman and Hall).

	  Also, have you asked about this directly to the maintainers of the 
"sn", "mnormt" and "mvtnorm" packages?  They might have other suggestions.

	  Hope this helps.
	  Spencer Graves
p.s.  Thanks for the self-contained example.  There seems to be a typo 
in your example:  Omega = diag(0, 27) is the matrix of all zeros, which 
produces a point mass at the center of the distribution.  I got your 
answers after changing it to diag(1, 27).

	  Making the dimension a variable, I found a sharp transition between k 
= 19 and 20:

 > k <- 19
 > xi <- alpha <- x <- rep(0,k)
 > Omega <- diag(1,k)
 > (p19 <- pmst(x, xi, Omega, alpha, df = 5))
[1] 1.907349e-06
attr(,"error")
[1] 2.413537e-20
attr(,"status")
[1] "normal completion"
 > k <- 20
 > xi <- alpha <- x <- rep(0,k)
 > Omega <- diag(1,k)
 > (p20 <- pmst(x, xi, Omega, alpha, df = 5))
[1] 0
attr(,"error")
[1] 1
attr(,"status")
[1] "oversize"

Konrad Banachewicz wrote:
> Dear All,
> I am using the pmst function from the sn package (version 0.4-0). After
> inserting the example from the help page, I get non-trivial answers, so
> everything is fine. However, when I try to extend it to higher dimension:
> xi <- alpha <- x <- rep(0,27)
> Omega <- diag(0,27)
> p1 <- pmst(x, xi, Omega, alpha, df = 5)
> 
> I get the following result:
> 
>> p1
> [1] 0
> attr(,"error")
> [1] 1
> attr(,"status")
> [1] "oversize"
> 
> So it seems like the dimension is a problem here (and not the syntax or type
> mismatch, as I inititally thought - the function is evaluated) - although I
> found no warning about it in the help page.
> 
> Can anyone give me a hint as to how to work around this problem and evaluate
> the skew-t cdf in a large-dimensional space? It's pretty crucial to my
> current research. Thanks in advance,
> 
> Konrad
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html