Simulation impliquant un conditionnement sur la somme de variables aléatoires

Je lisais cette question et j'ai pensé à simuler la quantité requise. Le problème est le suivant: si et sont normaux normaux, qu'est-ce que E (A ^ 2 | A + B) ? Je veux donc simuler E (A ^ 2 | A + B) . (pour une valeur choisie de A + B )$A$$B$$E(A^2|A+B)$$E(A^2|A+B)$$A+B$

J'ai essayé le code suivant pour y parvenir:

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 )

Le problème est qu'il n'y a presque toujours aucune valeur dans sum_ABlaquelle les correspondances x(entre les simulations). Si je choisis un élément sum_AB, c'est généralement la seule instance de sa valeur dans le vecteur.

En général, comment peut-on s'attaquer à ce problème et effectuer une simulation précise pour trouver une attente de la forme donnée? ( $A$ et $B$ ne sont pas nécessairement distribués normalement ou issus de la même distribution.)

Mon commentaire dans le fil référencé suggère une approche efficace: parce que et sont conjointement normaux avec une covariance nulle, ils sont indépendants, d'où la simulation n'a besoin que de générer (qui a une moyenne de et une variance ) et construire . Dans cet exemple, la distribution de est examinée au moyen de l'histogramme de valeurs simulées.$X=A+B$$Y=A-B$$Y$$0$$2$$A = (X+Y)/2$$A^2|(A+B=3)$$10^5$

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

L'attente peut être estimée comme

mean(a^2)

La réponse devrait être proche de .$11/4 = 2.75$

Une façon générique de résoudre ce problème consiste à considérer le changement de variables de à . Le jacobien de cette transformée étant égal à un (1), la densité de est Donc la densité de conditionnelle à est avec le terme de proportionnalité étant l'inverse de la densité marginale de , . Puisque , une transformée déterministe, c'est aussi la densité conjointe de donnée$(A,B)$$(A,A+B=S)$$(A,S)$

Vous pouvez résoudre ce problème en utilisant des exemples d'amorçage. Par exemple,

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)

En exécutant ce code par exemple, j'obtiens ce qui suit

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

Ainsi, lorsque alors .$A+B=0.9890429$$E(A^2|A+B=0.9890429)=0.7336328$

Maintenant pour valider que cela devrait être la réponse, exécutons le code de whuber dans sa solution. Donc, exécuter son code avec les x<-0.9890429résultats suivants:

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

Et donc les deux solutions sont très proches et coïncident l'une avec l'autre. Cependant, mon approche du problème devrait en fait vous permettre de saisir n'importe quelle distribution que vous souhaitez plutôt que de vous fier au fait que les données proviennent de distributions normales.

Une deuxième solution plus brutale qui repose sur le fait que lorsque la densité est relativement grande, vous pouvez facilement effectuer un calcul de force brute est la suivante

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

En exécutant ce code, nous obtenons ce qui suit

> estimate
[1] 2.757067

Ainsi, l'exécution du code pour donne ce qui correspond à la vraie solution.$A+B=3$$E(A^2|A+B=3)=2.757067$

il me semble que la question devient celle-ci:

1. comment simuler (X, Y) conditionnellement à X + Y = k puis
2. utiliser monte carlo pour estimer EU (X, Y) pour une fonction U (x, y)

Commençons par examiner l' échantillonnage d'importance :

$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)}$

où la première attente est relative à la variable aléatoire de densité et la seconde est par rapport à de densité .$Z_1$$f_1(z)$$Z_2$$f_2(z)$

Ainsi, si vous pouvez simuler aléatoirement les partir de puis estimer en utilisant ou alternativement simuler les partir de puis en utilisant$z_i$$f_1$$\frac{1}{n} \sum_i V(z_i)$$z_i$$f_2$$\frac{1}{n} \sum_i V(z_i) \frac{f_1(z_i)}{f_2(z_i)}$

Revenons maintenant à notre cas et sont distribués en tant que condition (X, Y) sur X + Y = k, c'est-à-dire et laissez$U(x,y)=x^2$$(X,Y)$$\frac{f(x,y)}{\int_{x+y=k} f(x,y)}$$A = \int_{x+y=k} f(x,y)$

alors maintenant la procédure est:

1. générer n copies iid à partir de la densité - et les appeler$g(x)$$X_i$
2. soit noter la distribution de ceci (X, Y) est , où est la fonction indicatrice$Y_i=k-X_i$$g(x)I(x+y=k)$$I()$
3. l'estimation est $$\frac{1}{n} \sum_i U(x_i,y_i) \frac{f(x_i,y_i)}{A g(x_i)}$$