[R] Fwd: dmnorm not meant for 1024-dimensional data?

Adelchi Azzalini aa at tango.stat.unipd.it
Wed Apr 25 21:49:59 CEST 2007

On Wed, Apr 25, 2007 at 01:25:09PM -0500, Daniel Elliott wrote:
> Hello.  Thank you for your mnormt package.  Below is an email I sent to
> R-help regarding the function dmnorm.
> Any help you can provide would be greatly appreciated.


You have not provided information on what specifically you 
have tried as for parameters and evaluation points, so I have
chosen a particularly simple case, namely the 1024-dimensional
density with independent standard marginals, using this code

d<- 1024
x <- 1
const<- 0.5*log(2*pi)
log.pdf <-0
for (i in 1:d)   log.pdf <- log.pdf - x^2/2 -const
S <- diag(d)
mu<- rep(0,d)
X <- rep(x,d)
log.pdf2 <- dmnorm(X,mu,S, log=TRUE)
cat(log.pdf, log.pdf2, abs(log.pdf-log.pdf2),"\n")

and the outcome was this one:
  -1453 -1453 2.137e-11 
which seems quite decent to me. Obviously, you must not take exp()
of this log-density, as otherwise a 0 is indeed produced, but 
the reason is not in the package, rather in the possibility of
representing that number using standard floating-point
numerical devices.

best regards,

Adelchi Azzalini

> Thank you.
> - dan elliott
> ---------- Forwarded message ----------
> From: Daniel Elliott <danelliottster at gmail.com>
> Date: Apr 25, 2007 1:21 PM
> Subject: dmnorm not meant for 1024-dimensional data?
> To: r-help at stat.math.ethz.ch
> Hello,
> I have some data generated by a simple mixture of Gaussians (more like
> K-means) and (as a test) am using dmnorm to calculate the probability
> of each data point coming from each Gaussian.  However, I get only
> zero probabilities.
> This code works in low dimensions (tried 2 and 3 already).  I have run
> into many implementations that do not work in high dimension, but I
> thought that I was safe with dmnorm because it has an option to
> compute the log of the probability.
> So, is dmnorm not intended to be used with data of such high dimensionality?
> Thank you,
> dan elliott

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