Es esto correcto ? (la generación de un Truncado-norma-multivariado de Gauss)

Question

Si $X\in\mathbb{R}^n,~X\sim \mathcal{N}(\underline{0},\sigma^2\mathbf{I})$ es decir, $$ f_X(x) = \frac{1}{{(2\pi\sigma^2)}^{n/2}} \exp\left(-\frac{||x||^2}{2\sigma^2}\right) $$

Quiero un análogo de la versión de un truncado-normal-distribución en un caso multivariante.

Más precisamente, quiero generar una norma-limitada a un valor de $\geq a$) multivariante de Gauss $Y$ s.t. $$ f_Y(y) = \begin{cases} c.f_X(y), \text{ if } ||y||\geq a \\[2mm] 0, \text{ otherwise }. \end{casos} $$ donde $c=\frac{1}{Prob\big\{||X||\geq a\big\}}$

---

Ahora observo lo siguiente:

Si $x=(x_1,x_2,\ldots,x_n)$, $||x||\geq a$

$\implies |x_n|\geq T\triangleq \sqrt{\max\left(0,\left(a^2-\sum_1^{n-1}x_i^2\right)\right)}$

Por lo tanto, por la elección de $x_1,\ldots,x_{n-1}$ como Gaussianas muestras, se puede restringir $x_n$ como una muestra de un Truncado-normal-distribución (después de una Gaussiana-cola $\geq T$) distribución $\mathcal{N}_T(0,\sigma^2)$, a excepción de su signo elegido aleatoriamente con una probabilidad de $1/2$.

Ahora mi pregunta es esta,

Si puedo generar cada vector de la muestra $(x_1,\ldots,x_n)$ $(X_1,\ldots,X_n)$ como

$x_1,\ldots,x_{n-1}\sim \mathcal{N}(0,\sigma^2)$

y

$x_n = Z_1 *Z_2~$ donde, $~Z_1\sim\{\pm1 ~\text{w.p.}~ 1/2\}$, $Z_2\sim\mathcal{N}_T(0,\sigma^2)$, (es decir, un truncado-escalar-normal RV con $T(x_1,\ldots,x_{n-1})\triangleq \sqrt{\max\left(0,\left(a^2-\sum_1^{n-1}x_i^2\right)\right)}$

Se $(X_1,X_2,\ldots,X_n)$ ser una norma-restringido ($\geq a$) multivariante de Gauss? (es decir, lo mismo que $Y$ definido anteriormente). ¿Cómo debo comprobar? Otras sugerencias si este no es el camino?

EDITAR:

Aquí es un scatter-plot de los puntos en 2D caso con la norma trunca a los valores por encima de "1"

Nota: Hay algunos grandes respuestas, pero la justificación de por qué esta propuesta es incorrecto falta. De hecho, eso es punto importante de esta pregunta.

Answer 1

La distribución normal multivariante de $X$ es esféricamente simétrica. La distribución que buscan trunca el radio de $\rho=||X||^2$ por debajo de lo $a$. Debido a que este criterio sólo depende de la longitud de $X$, el truncado distribución es esféricamente simétrica. Desde $\rho$ es independiente de la forma esférica ángulo de $X/||X||$ $\rho\,\sigma$ tiene un $\chi(n)$ distribución, por lo tanto, usted puede generar los valores de la distribución truncada en unos pocos y sencillos pasos:

1. Generar $X \sim \mathcal{N}(0,\mathbb{I}_n)$.

2. Generar $P$ como la raíz cuadrada de un $\chi^2(d)$ distribución truncada en $(a/\sigma)^2$.

3. Deje $Y = \sigma P\, X/||X||$.

En el paso 1, $X$ se obtiene como una secuencia de $d$ independiente realizaciones de una variable normal estándar.

En el paso 2, $P$ es fácilmente generados por la inversión de la función cuantil $F^{-1}$ $\chi^2(d)$ distribución: generar un uniforme de la variable $U$ apoyado en el rango (de cuantiles) entre $F((a/\sigma)^2)$ $1$ y establezca $P = \sqrt{F(U)}$.

Aquí es un histograma de $10^5$ independiente de la realización de $\sigma P$ $\sigma=3$ $n=11$ dimensiones, truncado por debajo de a $a=7$. Se tomó un segundo para generar, lo que demuestra la eficiencia del algoritmo.

La curva roja es la densidad de un truncado $\chi(11)$ distribución de la escala por $\sigma=3$. Su más cercano al histograma es la evidencia de la validez de esta técnica.

Para obtener una intuición para el truncamiento, considere el caso $a=3$, $\sigma=1$ en $n=2$ dimensiones. Aquí está un diagrama de dispersión de $Y_2$ contra $Y_1$ ($10^4$ independiente realizaciones). Esto demuestra claramente que el agujero en el radius $a$:

Por último, tenga en cuenta que (1) los componentes de $X_i$ debe tener distribuciones idénticas (debido a la simetría esférica) y (2) excepto cuando se $a=0$, que el común de distribución no es Normal. De hecho, como $a$ crece grande, la rápida disminución de la (univariante) distribución Normal, la causa de la mayoría de la probabilidad de que el esférico trunca multivariante normal clúster cerca de la superficie de la $n-1$-esfera (de radio $a$). La distribución marginal por lo tanto debe aproximarse a una escalada simétrica Beta$((n-1)/2,(n-1)/2)$ distribución concentrada en el intervalo de $(-a,a)$. Esto es evidente en el anterior diagrama de dispersión, donde $a=3\sigma$ ya es grande en dos dimensiones: los puntos de limón y un anillo (un $2-1$-esfera) de radio $3\sigma$.

Aquí están los histogramas de las distribuciones marginales a partir de una simulación de tamaño $10^5$ $3$ dimensiones con $a=10$, $\sigma=1$ (para que la aproximación de la Beta$(1,1)$ distribución es uniforme):

Desde el primer $n-1$ marginales del procedimiento descrito en la pregunta son normales (por construcción), que el procedimiento no puede ser correcta.

---

El siguiente R código generado en la primera figura. Se construye en paralelo de los pasos 1-3 para la generación de $Y$. Fue modificado para generar la segunda figura, cambiando variables a, d, ny sigma y la emisión de la trama de comando plot(y[1,], y[2,], pch=16, cex=1/2, col="#00000010") después y fue generado.

La generación de $U$ es modificado en el código de mayor resolución numérica: el código que realmente genera $1-U$ y la utiliza para calcular $P$.

La misma técnica de la simulación de los datos de acuerdo a una supuesta algoritmo, resumiendo, es con un histograma, y la superposición de un histograma puede ser utilizado para probar el método descrito en la pregunta. Va a confirmar que el método no funciona como se esperaba.

a <- 7      # Lower threshold
d <- 11     # Dimensions
n <- 1e5    # Sample size
sigma <- 3  # Original SD
#
# The algorithm.
#
set.seed(17)
u.max <- pchisq((a/sigma)^2, d, lower.tail=FALSE)
if (u.max == 0) stop("The threshold is too large.")
u <- runif(n, 0, u.max)
rho <- sigma * sqrt(qchisq(u, d, lower.tail=FALSE)) 
x <- matrix(rnorm(n*d, 0, 1), ncol=d)
y <- t(x * rho / apply(x, 1, function(y) sqrt(sum(y*y))))
#
# Draw histograms of the marginal distributions.
#
h <- function(z) {
  s <- sd(z)
  hist(z, freq=FALSE, ylim=c(0, 1/sqrt(2*pi*s^2)),
       main="Marginal Histogram",
       sub="Best Normal Fit Superimposed")
  curve(dnorm(x, mean(z), s), add=TRUE, lwd=2, col="Red")
}
par(mfrow=c(1, min(d, 4)))
invisible(apply(y, 1, h))
#
# Draw a nice histogram of the distances.
#
#plot(y[1,], y[2,], pch=16, cex=1/2, col="#00000010") # For figure 2
rho.max <- min(qchisq(1 - 0.001*pchisq(a/sigma, d, lower.tail=FALSE), d)*sigma, 
               max(rho), na.rm=TRUE)
k <- ceiling(rho.max/a)
hist(rho, freq=FALSE, xlim=c(0, rho.max),  
     breaks=seq(0, max(rho)+a, by=a/ceiling(50/k)))
#
# Superimpose the theoretical distribution.
#
dchi <- function(x, d) {
  exp((d-1)*log(x) + (1-d/2)*log(2) - x^2/2 - lgamma(d/2))
}
curve((x >= a)*dchi(x/sigma, d) / (1-pchisq((a/sigma)^2, d))/sigma, add=TRUE, 
      lwd=2, col="Red", n=257)

Answer 2

He escrito esto suponiendo que no desea que los puntos que tienen ||y||>, que es el análogo de la habitual en las dimensiones de truncamiento. Sin embargo, usted ha escrito que desea mantener los puntos que tienen |y|| >= a y desechamos los demás. Sin embargo, es obvio que el ajuste a mi de la solución puede ser hecho si usted realmente quiere mantener los puntos que tienen |y|| >= un.

La forma más sencilla, que pasa a ser un muy general técnica, es el uso de Aceptación-Rechazo https://en.wikipedia.org/wiki/Rejection_sampling . Va a ser bastante rápido mientras Prob(||X|| > a) es bastante baja, porque entonces no habrá muchos rechazos.

Generar una muestra valor de x a partir de las restricciones Multivariante Normal (aunque el problema de los estados que la Normal Multivariante es esférico, la técnica puede ser aplicada incluso si no lo está). Si ||x|| <= a aceptar, es decir, el uso de x, de lo contrario rechazarlo y generar una nueva muestra. Repita este proceso hasta que haya tantos aceptado muestras que usted necesita. El efecto de la aplicación de este procedimiento es generar y tal que su densidad es de c * f_X(y), si ||y|| <= a y 0 si ||y|| > una, por mi corrección a la apertura de parte de su pregunta. Usted nunca tendrá que calcular c; es, en efecto auto-determinada por el algoritmo basado en la frecuencia con la cual las muestras son rechazados.

Answer 3

Este es un buen intento pero no funciona debido a la "normalización constante": si consideramos el conjunto de la densidad $$f_X(x) \propto \frac{1}{{(2\pi\sigma^2)}^{n/2}} \exp\left(-\frac{||x||^2}{2\sigma^2}\right)\mathbb{I}_{||x||>a}=\frac{1}{{(2\pi\sigma^2)}^{n/2}} \exp\left(-\frac{x_1^2+\ldots+x_n^2}{2\sigma^2}\right)\mathbb{I}_{||x||>a}$$de la descomposición $$f_X(x) \propto \frac{1}{{(2\pi\sigma^2)}^{(n-1)/2}} \exp\left(-\frac{||x_{-n}||^2}{2\sigma^2}\right)\frac{1}{{(2\pi\sigma^2)}^{1/2}} \exp\left(-\frac{x_n^2}{2\sigma^2}\right)\mathbb{I}_{||x||>a}$$ $$=\frac{1}{{(2\pi\sigma^2)}^{(n-1)/2}} \exp\left(-\frac{||x_{-n}||^2}{2\sigma^2}\right)\frac{1}{{(2\pi\sigma^2)}^{1/2}} \exp\left(-\frac{x_n^2}{2\sigma^2}\right)\mathbb{I}_{||x_{-n}||^2+x_n^2>a^2}$$ $$=\frac{\mathbb{P}(X_n^2>a^2-||x_{-n}||^2)}{{(2\pi\sigma^2)}^{(n-1)/2}} \exp\left(-\frac{||x_{-n}||^2}{2\sigma^2}\right)\qquad\qquad\qquad\qquad\qquad$$ $$\qquad\qquad\qquad\times\frac{\mathbb{P}(X_n^2>a^2-||x_{-n}||^2)^{-1}}{{(2\pi\sigma^2)}^{1/2}} \exp\left(-\frac{x_n^2}{2\sigma^2}\right)\mathbb{I}_{x_n^2>a-||x_{-n}||^2}$$ que se integra a $$f_{X_{-n}}(x_{-n}) \propto \frac{\mathbb{P}(X_n^2>a^2-||x_{-n}||^2)}{{(2\pi\sigma^2)}^{(n-1)/2}} \exp\left(-\frac{||x_{-n}||^2}{2\sigma^2}\right)$$ en $x_n$, muestra que

1. La distribución condicional de $X_n$ dado los otros componentes, $X_{-n}$, es una distribución normal truncada;
2. La distribución marginal de los otros componentes, $X_{-n}$, es no una distribución normal debido a que el plazo adicional $\mathbb{P}(X_n^2>a^2-||x_{-n}||^2)$;

La única manera que puedo ver en tomar ventaja de esta propiedad es ejecutar un muestreador de Gibbs, uno de los componentes a la vez, el uso de la truncada normal distribuciones condicionales.

Answer 4

La pregunta se origina a partir de la idea de utilizar -- el condicionales básicas-descomposición de distribuciones conjuntas con el fin de dibujar el vector de muestras.

Deje $X$ ser un multivariante de Gauss con me.yo.d. componentes.

Deje $\text{Prob}(||X||>a) \triangleq T$ y $Y\triangleq X.\mathbb{I}_{||X||>a}$

El algoritmo en cuestión se propone sobre la base de los siguientes (todo correcto, pero engañando a la interpretación) condicional-factorización: $$f_Y(y) = \frac{1}{T}\frac{1}{{(2\pi\sigma^2)}^{n/2}} \exp\left(-\frac{||y||^2}{2\sigma^2}\right)\mathbb{I}_{||y||>a}\\ =\frac{1}{T}\frac{1}{{(2\pi\sigma^2)}^{n/2}} \exp\left(-\frac{y_1^2+\ldots+y_n^2}{2\sigma^2}\right)\mathbb{I}_{||y||>a}\\ =\left(\prod_{i=1}^{n-1} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{y_i^2}{2\sigma^2}\right)\right) \left(\frac{1}{T}\frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{y_n^2}{2\sigma^2}\right)\mathbb{I}_{||y||>a}\right)\\ =\underbrace{\left(\prod_{i=1}^{n-1} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{y_i^2}{2\sigma^2}\right)\right)}_{\text{Gaussianas}} \underbrace{\left(\frac{1}{T}\frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{y_n^2}{2\sigma^2}\right)\mathbb{I}_{y_n^2>(a^2-y_1^2-\ldots y_{n-1}^2)}\right)}_{\text{Gaussiano Truncado??}} $$

La respuesta más corta es que este último factor no es un truncado de Gauss, (lo más importante) ni siquiera una distribución.

---

Aquí está la explicación detallada de la razón por la que el de la factorización de la misma tiene un defecto fundamental. En una sola frase: condicional-factorización de un determinado conjunto de distribución deben satisfacer algunos de los fundamentales de las propiedades, y por encima de la factorización de no satisfacerlas (Ver más abajo).

En general, si alguna vez nos factorizar $f_{XY}(x,y)=f_X(x)\cdot f_{Y|X}(y|x)$ $f_X(x)$ es la marginal de $X$ $f_{Y|X}(y|x)$ es la distribución condicional de $Y$. Lo que significa:

1. El factor de $f(x,y)$ "se asume como" $f_X(x)$ debe ser de una distribución. Y,
2. El segundo factor "se asume como" $f_{Y|X}(y|x)$ debe ser una distribución para cada elección de $x$

En el ejemplo anterior, estamos tratando de condición como $Y_n|(Y_1\ldots Y_{n-1})$. Esto significa que la propiedad-1 para el factor de Gaussianas y propiedades-2 debe tener bueno para la última parte.

Es claro que la propiedad-1 se mantiene bien en el primer factor. Pero El problema es con la propiedad-2. El último factor anterior, desafortunadamente, no es una distribución (como el de olvidarse de Truncado de Gauss) para casi cualquier valor de $(Y_1\ldots Y_{n-1})$!!

---

Una propuesta de este tipo de algoritmo es probablemente un resultado de la siguiente equivocación: una Vez que una distribución natural de los factores de una distribución conjunta (como Gaussianas en la anterior), que conduce a una condicional de la factorización. \begin{cases} c.f_X(y), \text{ if } ||y||\geq a \\[2mm] 0, \text{ otherwise }. \endNo! ---- El otro (segundo) factor también debe ser bueno.

---

Nota: No es una gran respuesta por parte de @whuber anteriores, que en realidad no resuelve el problema de la generación de una norma trunca multivariante de Gauss. Estoy aceptando su respuesta. Esta respuesta es sólo para aclarar y compartir mi propio entendimiento y la génesis de la cuestión.