# Simulation involving conditioning on sum of random variables

I was reading this question, and thought about simulating the required quantity. The problem is as follows: If $A$ and $B$ are iid standard normal, what is $E(A^2|A+B)$? So I want to simulate $E(A^2|A+B)$. (for a chosen value of $A+B$)

I tried the following code to achieve this:

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 )


The problem is that there is almost always no value in sum_AB which matches x (across simulations). If I choose some element from sum_AB, then it usually the only instance of its value in the vector.

In general, how can one tackle this problem and perform an accurate simulation to find an expectation of the given form? ($A$ and $B$ may not necessarily be normally distributed, or from the same distribution.)

-
Your recent edit substantially changes the question, as our interchange of comments indicates. It becomes more difficult to answer in the much greater generality you now suppose. For instance, there exist special--and rather involved--techniques just to answer it when the value of $A+B$ is rare (out in one of the tails). –  whuber Aug 25 at 21:24
@whuber Wouldn't all values be relatively rare when we are dealing with two continuous random variables? –  Comp_Warrior Aug 25 at 21:30
Yes, but narrow bands of values--which usually suffice for such simulations--would never work out in the tails (nor in any other region where the PDF gets very small), whereas when the density is relatively large you can easily perform a brute-force calculation that is assured of producing a decent number of data having $A+B$ close enough to its desired value to enable some conclusions to be drawn from the simulation. –  whuber Aug 25 at 21:32
@whuber I see - could you give some indication in your answer of the special techniques you mention? Apologies for not indicating what I was interested in below in the comments. –  Comp_Warrior Aug 25 at 21:34
Comp_Warrior I am appending a second solution which I believe is what @whuber is alluding to. –  Dan Aug 25 at 21:48

My comment in the referenced thread suggests one efficient approach: because $X=A+B$ and $Y=A-B$ are jointly Normal with zero covariance, they are independent, whence the simulation only needs to generate $Y$ (which has mean $0$ and variance $2$) and construct $A = (X+Y)/2$. In this example the distribution of $A^2|(A+B=3)$ is examined by means of the histogram of $10^5$ simulated values.

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


The expectation can be estimated as

mean(a^2)


The answer should be close to $11/4 = 2.75$.

-
Thanks - this makes sense. However, am I right in understanding that this simplification will only work if both random variables in question are iid normal? What if I had a case where $A$ and $B$ were from another distribution (and possible separate from each other)? –  Comp_Warrior Aug 25 at 20:56
Your understanding is correct. This is one reason Normal variables are so popular, both theoretically and in computer models! Nevertheless, the basic idea of looking for a way to transform the variables into sets of independent (or easily related) variables will carry over to a more general setting. –  whuber Aug 25 at 20:57

You could solve this problem using bootstrap samples. For example,

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)


Running this code for example, I obtain the following

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


So when $A+B=0.9890429$ then $E(A^2|A+B=0.9890429)=0.7336328$.

Now to validate that this should be the answer, let's run whuber's code in his solution. So running his code with x<-0.9890429 results in the following:

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


And so the two solutions are very close and coincide with one another. However, my approach to the problem should actually allow you to input any distribution you want rather than relying on the fact that the data came from Normal distributions.

A second more so brute force solution that relies on the fact that when the density is relatively large you can easily perform a brute-force calculation is the following

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


Running this code we obtain the following

> estimate
[1] 2.757067


Thus running the code for $A+B=3$ results in $E(A^2|A+B=3)=2.757067$ which agrees with the true solution.

-
I must be missing something: the question asks for the user to specify the value of $A+B$. Where is that done in your code? What would your code look like in the case $A+B$ needed to be set to $3$, for instance? –  whuber Aug 25 at 21:22
@whuber you are completely correct. I can only do it for sums that I know will appear. –  Dan Aug 25 at 21:24

it seems to me that the question becomes this:

1. how to simulate (X,Y) conditional on X+Y=k and then
2. use monte carlo to estimate EU(X,Y) for some function U(x,y)

let's start by reviewing importance sampling :

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

where the first expectations is with respect to random variable $Z_1$ with density $f_1(z)$ and the second one is wrt $Z_2$ with density $f_2(z)$.

Thus if you can randomly simulate $z_i$'s from $f_1$ then estimate using $\frac{1}{n} \sum_i V(z_i)$ or alternatively simulate $z_i$'s from $f_2$ then using $\frac{1}{n} \sum_i V(z_i) \frac{f_1(z_i)}{f_2(z_i)}$

Now let's get back to our case $U(x,y)=x^2$ and $(X,Y)$ are distributed as (X,Y) condition on X+Y=k, i.e. $\frac{f(x,y)}{\int_{x+y=k} f(x,y)}$ and let $A = \int_{x+y=k} f(x,y)$

so now the procedure is :

1. generate n iid copies from density $g(x)$ - and call them $X_i$
2. let $Y_i=k-X_i$ note the distribution of this (X,Y) is $g(x)I(x+y=k)$, where $I()$ is indicator function
3. the estimate is $$\frac{1}{n} \sum_i U(x_i,y_i) \frac{f(x_i,y_i)}{A g(x_i)}$$
-