La simulación involucra acondicionado en la suma de variables aleatorias

Question

Estaba leyendo esta pregunta, y el pensamiento acerca de la simulación de la cantidad requerida. El problema es el siguiente: Si $A$ $B$ son iid normal estándar, ¿cuál es $E(A^2|A+B)$? Así que quiero simular $E(A^2|A+B)$. (para un valor seleccionado de $A+B$)

He probado el siguiente código para lograr esto:

n <- 1000000
x <- 1 # the sum of A and B

A <- rnorm(n)
B <- rnorm(n)

sum_AB = A+B

estimate <- 1/sum(sum_AB==x) * sum( (A[sum_AB==x])^2 )

El problema es que casi siempre hay ningún valor en sum_AB que coincide con x (a través de simulaciones). Si he de elegir algún elemento de sum_AB, entonces es generalmente la única instancia de su valor en el vector.

En general, ¿cómo se puede abordar este problema y realizar una simulación precisa para encontrar una expectativa de la forma determinada? ($A$$B$ puede no ser necesariamente una distribución normal, o de la misma distribución.)

Answer 1

Mi comentario en el que se hace referencia hilo sugiere un enfoque eficaz: porque $X=A+B$ $Y=A-B$ son conjuntamente Normal con cero covarianza, son independientes, donde la simulación sólo necesita generar $Y$ (que ha decir $0$ y la varianza $2$) y construcción $A = (X+Y)/2$. En este ejemplo, la distribución de $A^2|(A+B=3)$ es examinado por medio de un histograma de $10^5$ valores simulados.

x <- 3
y <- rnorm(1e5, 0, sqrt(2))
a <- (x+y)/2
hist(a^2)

La expectativa puede ser estimado como

mean(a^2)

La respuesta debe estar cerca de la $11/4 = 2.75$.

Answer 2

Usted puede resolver este problema usando bootstrap muestras. Por ejemplo,

n <- 1000000

A <- rnorm(n)
B <- rnorm(n)
AB <- cbind(A,B)

boots <- 100
bootstrap_data <- matrix(NA,nrow=boots*n,ncol=2)


for(i in 1:boots){
    index <- sample(1:n,n,replace=TRUE)
    bootstrap_data[(i*n-n+1):(i*n),] <- cbind(A[index],B[index]) 
}

sum_AB <- bootstrap_data[,1] + bootstrap_data[,2]
x <- sum_AB[sample(1:n,1)]

idx <- which(sum_AB == x)

estimate <- mean(bootstrap_data[idx,1]^2)

Ejecutando este código, por ejemplo, obtener el siguiente

> estimate
[1] 0.7336328
> x
[1] 0.9890429

Así que cuando $A+B=0.9890429$$E(A^2|A+B=0.9890429)=0.7336328$.

Ahora para validar que esta debe ser la respuesta, vamos a correr whuber del código en su solución. Por lo que se ejecuta su código con x<-0.9890429 de resultados en las siguientes:

> x <- 0.9890429
> y <- rnorm(1e5, 0, sqrt(2))
> a <- (x+y)/2
> hist(a^2)
>
> mean(a^2)
[1] 0.745045

Y así las dos soluciones están muy cerca y coinciden el uno con el otro. Sin embargo, mi enfoque para el problema en realidad debería permitir a la entrada de cualquier distribución que usted desea en lugar de confiar en el hecho de que los datos provenían de distribuciones Normales.

---

Un segundo más para solución de fuerza bruta que se basa en el hecho de que cuando la densidad es relativamente grande, usted puede realizar un ataque de fuerza bruta de cálculo es la siguiente

n <- 1000000

x <- 3  #The desired sum to condition on

A <- rnorm(n)
B <- rnorm(n)
sum_AB <- A+B

epsilon <- .01
idx <- which(sum_AB > x-epsilon & sum_AB < x+epsilon)
estimate <- mean(A[idx]^2)

estimate

Ejecutando este código, obtenemos los siguientes

> estimate
[1] 2.757067

Por lo tanto se ejecuta el código para $A+B=3$ resultados en $E(A^2|A+B=3)=2.757067$ que está de acuerdo con la verdadera solución.

Answer 3

a mí me parece que la pregunta es esta:

1. cómo simular (X,Y) condicional en X+Y=k y, a continuación,
2. el uso de monte carlo para la estimación de la UE(X,Y) para alguna función U(x,y)

vamos a empezar por la revisión de la importancia de muestreo :

$E V(Z_1) = \int V(z) f_1(z) = \int V(z) \frac{f_1(z)}{f_2(z)} f_2(z) = E V(Z_2)\frac{f_1(Z_2)}{f_2(Z_2)}$

donde el primer expectativas con respecto a la variable aleatoria $Z_1$ con densidad de $f_1(z)$ y el segundo es wrt $Z_2$ con densidad de $f_2(z)$.

Por lo tanto, si al azar se puede simular $z_i$'s de $f_1$ luego estimar el uso de $\frac{1}{n} \sum_i V(z_i)$ o, alternativamente, simular $z_i$'s de $f_2$ uso de $\frac{1}{n} \sum_i V(z_i) \frac{f_1(z_i)}{f_2(z_i)}$

Ahora volvamos a nuestro caso $U(x,y)=x^2$ $(X,Y)$ se distribuye de la forma (X,Y) en la condición X+Y=k, es decir, $\frac{f(x,y)}{\int_{x+y=k} f(x,y)}$ y deje $A = \int_{x+y=k} f(x,y)$

así que ahora el procedimiento es el siguiente :

1. generar n iid copias de densidad de $g(x)$ - y llamarlos $X_i$
2. deje $Y_i=k-X_i$ nota: la distribución de esta (X,Y) es $g(x)I(x+y=k)$ donde $I()$ es indicador de la función
3. la estimación es $$\frac{1}{n} \sum_i U(x_i,y_i) \frac{f(x_i,y_i)}{A g(x_i)} $$