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:
Generar $X \sim \mathcal{N}(0,\mathbb{I}_n)$.
Generar $P$ como la raíz cuadrada de un $\chi^2(d)$ distribución truncada en $(a/\sigma)^2$.
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.
![Figure]()
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$:
![Figure 2]()
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):
![Figure 3]()
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
, n
y 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)