Generar muestras de singular distribución Gausiana

Question

Vamos aleatorios vectoriales $x = (x_1,...,x_n)$ seguir multivariante distribución normal con media de $m$ y matriz de covarianza $S$. Si $S$ es simétrica y definida positiva (que es el caso habitual), entonces uno puede generar muestras aleatorias de $x$ por el primer muestreo indepently $r_1,...,r_n$ de la normal estándar y, a continuación, utilizando la fórmula $m + Lr$ donde $L$ es el Cholesky menor factor para $S=LL^T$$r = (r_1,...,r_n)^T$.

¿Qué pasa si uno quiere que las muestras de singular Gaussiano es decir $S$ todavía es simétrica pero no más positiva definida (sólo positiva semi-definido). Podemos suponer también que las desviaciones (diagonal elementos de $S$) son estrictamente positivos. A continuación, algunos elementos de la $x$ debe tener relación lineal y la distribución en realidad se encuentra en menor espacio tridimensional con la dimensión de $<n$, ¿verdad?

Es obvio que si, por ejemplo,$n=2, m = \begin{bmatrix} 0 \\ 0 \end{bmatrix}, S = \begin{bmatrix} 1 & 1 \\ 1 & 1\end{bmatrix}$, entonces uno puede generar $x_1 \sim N(0,1)$ y establezca $x_2=x_1$ ya que están totalmente correlacionados. Sin embargo, hay una buena métodos para la generación de muestras para el caso general,$n>2$? Supongo que uno debe ser capaz de identificar las dimensiones inferiores subespacio, a continuación, pasar a ese espacio en donde uno será válido matriz de covarianza, a continuación, muestra de ello y, finalmente, deducir los valores de la linealmente dependiente de las variables de este menor dimensiones de la muestra. Pero ¿cuál es la mejor manera para que en la práctica? Puede que alguien me señale a los libros o artículos que se ocupan del tema; yo no podía encontrar uno.

Answer 1

La singular distribución de Gauss, es el impulso hacia adelante de un nonsingular de distribución en un menor espacio tridimensional. Geométricamente, usted puede tomar una distribución Normal estándar, escalarlo, rotarlo, y se incrusta isométricamente en una afín subespacio de un espacio de dimensiones superiores. Algebraicamente, esto se hace por medio de una Descomposición de Valor Singular (SVD) o su equivalente.

---

Deje $\Sigma$ ser la matriz de covarianza y $\mu$ la media en $\mathbb{R}^n$. Debido a $\Sigma$ es no negativa definida y simétrica, la enfermedad vesicular porcina se toman la forma

$$\Sigma = U \Lambda^2 U^\prime$$

for an orthogonal matrix $U\en O(n)$ and a diagonal matrix $\Lambda$. $\Lambda$ will have $m$ nonzero entries, $0\le m \le n$.

Let $X$ have a standard Normal distribution in $\mathbb{R}^m$: that is, each of its $m$ components is a standard Normal distribution with zero mean and unit variance. Abusing notation a little, extend the components of $X$ with $n-m$ zeros to make it an $n$-vector. Then $U\Lambda X$ is in $\mathbb{R}^n$ and we may compute

$$\text{Cov}(U\Lambda X) = U \Lambda\text{Cov}(X) \Lambda^\prime U^\prime = U \Lambda^2 U^\prime = \Sigma.$$

Consequently

$$Y = \mu + U\Lambda X$$

has the intended Gaussian distribution in $\mathbb{R}^n$.

It is of interest that this works when $n=m$: that is to say, this is a (standard) method to generate multivariate Normal vectors, in any dimension, for any given mean $\mu$ and covariance $\Sigma$ by using a univariate generator of standard Normal values.

---

As an example, here are two views of a thousand simulated points for which $n=3$ and $m=2$:

El segundo punto de vista, a partir de la orilla -, demuestra la singularidad de la distribución. El R código que producen estas cifras sigue a la anterior matemático de la exposición.

#
# Specify a Normal distribution.
#
mu <- c(5, 5, 5)
Sigma <- matrix(c(1, 2, 1,
                  2, 3, 1,
                  1, 1, 0), 3)
#
# Analyze the covariance.
#
n <- dim(Sigma)[1]
s <- svd((Sigma + t(Sigma))/2) # Guarantee symmetry
s$d <- abs(zapsmall(s$d))
m <- sum(s$d > 0)
    #$
# Generate a standard Normal `x` in R^m.
#
n.sample <- 1e3 # Number of points to generate
x <- matrix(rnorm(m*n.sample), nrow=m)
#
# Embed `x` in R^n and apply the square root of Sigma obtained from its SVD.
#
x <- rbind(x, matrix(0, nrow=n-m, ncol=n.sample))
y <- s$u %*% diag(sqrt(s$d)) %*% x + mu
#
# Plot the results (presuming n==3).
#
library(rgl)
plot3d(t(y), type="s", size=1, aspect=TRUE, 
       xlab="Y1", ylab="Y2", zlab="Y3", box=FALSE,
       col="Orange")