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