¿Cómo calcular la probabilidad total dentro de un corte de una distribución normal bivariada en R?

Question

Tengo una distribución normal bivariada compuesta por las distribuciones normales univariadas $X_1$ y $X_2$ con $\rho \approx 0.3$ .

$$ \begin{pmatrix} X_1 \\ X_2 \end{pmatrix} \sim \mathcal{N} \left( \begin{pmatrix} \mu_1 \\ \mu_2 \end{pmatrix} , \begin{pmatrix} \sigma^2_1 & \rho \sigma_1 \sigma_2 \\ \rho \sigma_1 \sigma_2 & \sigma^2_2 \end{pmatrix} \right) $$

¿Existe una forma sencilla de calcular en R la probabilidad acumulada de $X_1$ siendo menor que un valor $z$ dada una porción particular de $X_2$ (entre dos valores $a,b$ ) dado que conocemos todos los parámetros $\mu_1, \mu_2, \sigma_1, \sigma_2, \rho$ ?

$P(X_1 < z | a < X_2 < b)$

¿Puede la función de distribución que busco coincidir (o ser aproximada por) la función de distribución de una distribución normal univariante (para usar qnorm / pnorm )? Lo ideal sería que fuera así para poder realizar el cálculo con menos dependencias de bibliotecas (por ejemplo, en un servidor MySQL).

Esta es la distribución bivariada que estoy utilizando:

means <- c(79.55920, 52.29355)
variances <- c(268.8986, 770.0212)
rho <- 0.2821711

covariancePartOfMatrix <- sqrt(variances[1]) * sqrt(variances[2]) * rho
sigmaMatrix <- matrix(c(variances[1],covariancePartOfMatrix,covariancePartOfMatrix,variances[2]), byrow=T, ncol=2)

n <- 10000
dat <- MASS::mvrnorm(n=n, mu=means, Sigma=sigmaMatrix)

plot(dat)

Este es mi intento numérico de obtener el resultado correcto. Sin embargo, utiliza datos generados de la distribución bivariada y no estoy convencido de que dé el resultado correcto.

a <- 79.5
b <- 80.5
z <- 50

sliceOfDat <- subset(data.frame(dat), X1 > a, X1 < b)
estimatedMean <- mean(sliceOfDat[,c(2)])
estimatedDev <- sd(sliceOfDat[,c(2)])

estimatedPercentile <- pnorm(z, estimatedMean, estimatedDev)

Edición - Implementación en R de la solución basada en la respuesta de whuber

A continuación se presenta una implementación de la solución aceptada utilizando integrate comparado con mi idea original basada en el muestreo. La solución aceptada proporciona el resultado esperado 0,5, mientras que mi idea original se desviaba en una cantidad significativa (0,41). Actualización - Ver la edición de wheber para una mejor implementación.

# Bivariate distribution parameters
means <- c(79.55920, 52.29355)
variances <- c(268.8986, 770.0212)
rho <- 0.2821711

# Generate sample data for bivariate distribution
n <- 10000

covariancePartOfMatrix <- sqrt(variances[1]) * sqrt(variances[2]) * rho
sigmaMatrix <- matrix(c(variances[1],covariancePartOfMatrix,covariancePartOfMatrix,variances[2]), byrow=T, ncol=2)
dat <- MASS::mvrnorm(n=n, mu=means, Sigma=sigmaMatrix)

# Input parameters to test the estimation
w = 79.55920

a <- w - 0.5
b <- w + 0.5
z <- 52.29355

# Univariate approximation using randomness
sliceOfDat <- subset(data.frame(dat), X1 > a, X1 < b)
estimatedMean <- mean(sliceOfDat[,c(2)])
estimatedDev <- sd(sliceOfDat[,c(2)])

estimatedPercentile <- pnorm(z, estimatedMean, estimatedDev)
# OUTPUT: 0.411

# Numerical approximation from exact solution
adaptedZ <- (z - means[2]) / sqrt(variances[2])
adaptedB <- (b - means[1]) / sqrt(variances[1])
adaptedA <- (a - means[1]) / sqrt(variances[1])

exactSolutionCoeff <- 1 / (pnorm(adaptedB) - pnorm(adaptedA))
integrand <- function(x) pnorm((adaptedZ - rho * x) / sqrt(1 - rho * rho)) * dnorm(x)
exactSolutionInteg <- integrate(integrand, adaptedA, adaptedB)
# 0.0121, abs.error 1.348036e-16, "OK"
exactPercentile = exactSolutionCoeff * exactSolutionInteg$value
# OUTPUT: 0.500

Answer 1

Sí, una aproximación normal funciona, pero no en todos los casos. Tenemos que hacer algunos análisis para identificar cuándo la aproximación es buena.

Solución exacta

Reexpresar $(X_1,X_2)$ en unidades estandarizadas para que tengan medias y varianzas unitarias nulas. Dejando que $\Phi$ sea la función de distribución normal estándar (su FCD), es bien sabido por la teoría de la regresión por mínimos cuadrados ordinarios que

$$\Pr(X_1 \le z\,|\, X_2 = x) = \Phi\left(\frac{z - \rho x}{\sqrt{1-\rho^2}}\right).$$

La probabilidad deseada puede obtenerse entonces mediante la integración:

$$\eqalign{\Pr(X_1 \le z\,|\, a \lt X_2 \le b) &= \frac{1}{\Phi(b)-\Phi(a)}\int_a^b \Pr(X_1\le z\,|\, X_2=x) \phi(x)\,dx \\&= \frac{1}{\Phi(b)-\Phi(a)}\int_a^b \Phi\left(\frac{z - \rho x}{\sqrt{1-\rho^2}}\right) \phi(x)\,dx.}$$

Esto parece requerir una integración numérica (aunque el resultado para $(a,b)=\mathbb{R}$ se puede obtener de forma cerrada: véase ¿Cómo puedo calcular $\int^{\infty}_{-\infty}\Phi\left(\frac{w-a}{b}\right)\phi(w)\,\mathrm dw$ ).

Aproximación de la distribución

Esta expresión se puede diferenciar bajo el signo integral (con respecto a $z$ ) para obtener el PDF,

$$f(z\,|\, a \lt X_2 \le b) = \phi(z)\ \frac{\Phi\left(\frac{b-\rho z}{\sqrt{1-\rho^2}}\right) - \Phi\left(\frac{a-\rho z}{\sqrt{1-\rho^2}}\right)}{\Phi(b) - \Phi(a)}.$$

Esto muestra el PDF como producto del PDF normal estándar $\phi$ y una "corrección". Cuando $b-a$ es pequeño en comparación con $\sqrt{1-\rho^2}$ (concretamente, cuando $(b-a)^2 \ll 1-\rho^2$ ), podríamos aproximar la diferencia en el numerador con la primera derivada:

$$\Phi\left(\frac{b-\rho z}{\sqrt{1-\rho^2}}\right) - \Phi\left(\frac{a-\rho z}{\sqrt{1-\rho^2}}\right)\approx \phi\left(\frac{(a+b)/2-\rho z}{\sqrt{1-\rho^2}}\right)\frac{b-a}{\sqrt{1-\rho^2}}.$$

El error de esta aproximación está uniformemente acotado (en todos los valores de $z$ ) porque la segunda derivada de $\Phi$ está acotado.

Con esta aproximación, y completando el cuadrado, obtenemos

$$f(z\,|\, a \lt X_2 \le b) \approx \phi\left(z; \rho(a+b)/2, \sqrt{1-\rho^2}\right) \frac{(b-a)\exp\left(-(a+b)^2/8\right)}{(\Phi(b)-\Phi(a))\sqrt{2\pi}}.$$

( $\phi(*; \mu, \sigma)$ denota la FDP de una distribución normal de media $\mu$ y la desviación estándar $\sigma$ .)

Además, el factor de la derecha (que no depende de $z$ ) debe estar muy cerca de $1$ porque el $\phi$ se integra a la unidad, por lo que

$$f(z\,|\, a \lt X_2 \le b) \approx \phi\left(z; \rho(a+b)/2, \sqrt{1-\rho^2}\right).$$

Todo esto tiene mucho sentido: la distribución condicional de $X_1$ se aproxima por su distribución condicional en el punto medio del intervalo, $(a+b)/2$ , donde tiene media $\rho(a+b)/2$ y la desviación estándar $\sqrt{1-\rho^2}$ . El error es proporcional a la anchura del intervalo $b-a$ y a una expresión dominada por $\exp(-(a+b)^2/8)$ que sólo es importante cuando ambos $a$ y $b$ están fuera en la misma cola. Por lo tanto, esta aproximación normal funciona para rodajas estrechas no demasiado lejos en las colas de la distribución bivariada. Además, la diferencia entre

$$\frac{(b-a)\exp\left(-(a+b)^2/8\right)}{(\Phi(b)-\Phi(a))\sqrt{2\pi}}$$

y $1$ sirve para comprobar la calidad de la aproximación.

---

Editar

Para comprobar estas conclusiones, he simulado datos en R para varios valores de $b$ y $\rho$ ( $a=-3$ en todos los casos), dibujé su densidad empírica y superpuse sobre ella las densidades teóricas (azul) y aproximadas (rojo) para compararlas. (No se pueden ver los gráficos de densidad porque los gráficos teóricos se ajustan a ellos casi perfectamente). Como $|\rho|$ se acerca a $1$ la aproximación se vuelve más pobre: esto merece un estudio más profundo. Claro que la aproximación es excelente para valores suficientemente pequeños de $b-a$ .

#
# Numerical integration, to give a correct value.
#
f <- function(z, a, b, rho, value.only=FALSE, ...) {
  g <- function(x) pnorm((z - rho*x)/sqrt(1-rho^2)) * dnorm(x) / (pnorm(b) - pnorm(a))
  u <- integrate(g, a, b, ...)
  if (value.only) return(u$value) else return(u)
}
#
# Set up the problem.
#
a <- -3  # Left endpoint
n <- 1e5 # Simulation size
par(mfrow=c(2,3))
for (rho in c(1/4, -2/3)) {
  for (b in c(-2.5, -2, -1.5)) {
    z <- seq((a-3)*rho, (b+3)*rho, length.out=101)
    #
    # Check the approximation (`v` needs to be small).
    #
    v <- (b-a) * exp(-(a+b)^2/8) / (pnorm(b) - pnorm(a)) / sqrt(2*pi) - 1
    #
    # Simulate some values of (x1, x2).
    #
    x.2 <- qnorm(runif(n, pnorm(a), pnorm(b)))  # Normal between `a` and `b`
    x.1 <- rho*x.2 + rnorm(n, sd=sqrt(1-rho^2))
    #
    # Compare the simulated distribution to the theoretical and approximate
    # densities.
    #
    x.hat <- density(x.1)
    plot(x.hat, type="l", lwd=2, 
         main="Simulated, True, and Approximate",
         sub=paste0("a=", round(a,2), ", b=", round(b, 2), ", rho=", round(rho, 2), 
         "; v=", round(v,3)),
         xlab="X.1", ylab="Density")

    # Theoretical
    curve(dnorm(x) * (pnorm((b-rho*x)/sqrt(1-rho^2)) - pnorm((a-rho*x)/sqrt(1-rho^2))) /
       (pnorm(b) - pnorm(a)), lwd=2, col="Blue", add=TRUE)

    # Approximate
    curve(dnorm(x, rho*(a+b)/2, sqrt(1-rho^2)), col="Red", lwd=2, add=TRUE)
  }
}