[R] Multivariate skew-t cdf
Adelchi Azzalini
aa at tango.stat.unipd.it
Mon Jun 5 19:36:01 CEST 2006
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/
More information about the R-help
mailing list