Simulación que involucra condicionamiento sobre la suma de variables aleatorias

Estaba leyendo esta pregunta y pensé en simular la cantidad requerida. El problema es el siguiente: si y son estándar normal, ¿qué es ? Entonces quiero simular . (para un valor elegido de )$A$$B$$E(A^2|A+B)$$E(A^2|A+B)$$A+B$

Intenté 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 no hay ningún valor en el sum_ABque coincidan x(entre simulaciones). Si elijo algún elemento sum_AB, generalmente es 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 dada? ( y no necesariamente se distribuyen normalmente, o de la misma distribución).$A$$B$

Mi comentario en el hilo de referencia sugiere un enfoque eficiente: porque e son conjuntamente normales con cero covarianza, son independientes, por lo que la simulación solo necesita generar (que tiene una media de$X=A+B$$Y=A-B$$Y$$0$ y varianza$2$) y construir $A = (X+Y)/2$. En este ejemplo, la distribución de$A^2|(A+B=3)$ se examina mediante el histograma de $10^5$ valores simulados

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

La expectativa se puede estimar como

mean(a^2)

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

Una forma genérica de resolver este problema es considerar el cambio de variables de $(A,B)$ a $(A,A+B=S)$. El jacobiano de esta transformación es igual a uno (1), la densidad de$(A,S)$ es

$$f_{A,S}(a,s)=f_A(a)f_B(s-a)$$ Por lo tanto, la densidad de $A$ condicional en $S=s$ es $$f_{A|S}(a|s)\propto f_A(a)f_B(s-a)$$ siendo el término de proporcionalidad el inverso de la densidad marginal de $S$, $f_S(s)^{-1}$. Ya que$B=S-A$, una transformación determinista, esta es también la densidad conjunta de $(A,B)$ dado $S$ $$f_{A,B|S}(a,b|s)\propto f_A(a)f_B(s-a)\mathbb{I}_{a+b=s}$$ La generación de una realización a partir de este objetivo se puede hacer directamente si la forma es lo suficientemente simple, o por aceptar-rechazar, Metropolis-Hastings, muestreo de corte o cualquier otro método de simulación estándar.

Puede resolver este problema utilizando muestras de bootstrap. 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)

Al ejecutar este código, por ejemplo, obtengo lo siguiente

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

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

Ahora para validar que esta debería ser la respuesta, ejecutemos el código de Whuber en su solución. Así que ejecuta su código con x<-0.9890429resultados en lo siguiente:

> 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 son muy cercanas y coinciden entre sí. Sin embargo, mi enfoque del problema debería permitirle ingresar cualquier distribución que desee en lugar de confiar en el hecho de que los datos provienen de distribuciones normales.

---

Una segunda solución más de fuerza bruta que se basa en el hecho de que cuando la densidad es relativamente grande, puede realizar fácilmente un cálculo de fuerza bruta 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

Al ejecutar este código obtenemos lo siguiente

> estimate
[1] 2.757067

Ejecutando así el código para $A+B=3$ resultados en $E(A^2|A+B=3)=2.757067$ que concuerda con la verdadera solución.

me parece que la pregunta se convierte en esto:

1. cómo simular (X, Y) condicional en X + Y = k y luego
2. use monte carlo para estimar EU (X, Y) para alguna función U (x, y)

Comencemos por revisar el muestreo de importancia :

$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 las primeras expectativas son con respecto a la variable aleatoria $Z_1$ con densidad $f_1(z)$ y el segundo es wrt $Z_2$ con densidad $f_2(z)$.

Por lo tanto, si puede simular aleatoriamente $z_i$es de $f_1$ luego estimar usando $\frac{1}{n} \sum_i V(z_i)$ o alternativamente simular $z_i$es de $f_2$ luego usando $\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$ y $(X,Y)$ se distribuyen como (X, Y) condición en X + Y = k, es decir $\frac{f(x,y)}{\int_{x+y=k} f(x,y)}$ y deja $A = \int_{x+y=k} f(x,y)$

así que ahora el procedimiento es:

1. generar copias n iid de densidad $g(x)$ - y llamarlos $X_i$
2. dejar $Y_i=k-X_i$ tenga en cuenta que la distribución de este (X, Y) es $g(x)I(x+y=k)$, dónde $I()$ es la función del indicador
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)}$$