[R] Multivariate skew-t cdf

[Adelchi Azzalini]

On Mon, Jun 05, 2006 at 09:43:24AM -0700, Spencer Graves wrote:


Thanks to Spencer Graves for provinding this clarification about pmst.
As the author of pmst, I was really the one expected to answer the
query, but  I was on travel in the last two weeks and did not read 
this query.

As a complement to what already explained, the reason of the sharp change
from k=19 to k=20 is as follows: pmst (package sn) calls pmt (package
mnormt) with k increased by 1, and pmt calls a Fortran routine (written by
Alan Genz), which issues an error if k>20 or k<1. Hence, the effective
maximum number of dimensions allowed by pmst is 19. 

best wishes,

Adelchi Azzalini



> 	  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
> 
> ______________________________________________
> 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
> 

-- 
Adelchi Azzalini  <azzalini at stat.unipd.it>
Dipart.Scienze Statistiche, Università di Padova, Italia
tel. +39 049 8274147,  http://azzalini.stat.unipd.it/