Simulation involving conditioning on sum of random variables

Question

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.)

Answer 1

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$.

Answer 2

A generic way to solve this problem is to consider the change of variables from $(A,B)$ to $(A,A+B=S)$. The Jacobian of this transform being equal to one (1), the density of $(A,S)$ is $$f_{A,S}(a,s)=f_A(a)f_B(s-a)$$ Therefore the density of $A$ conditional on $S=s$ is $$f_{A|S}(a|s)\propto f_A(a)f_B(s-a)$$ with the proportionality term being the inverse of the marginal density of $S$, $f_S(s)^{-1}$. Since $B=S-A$, a deterministic transform, this is also the joint density of $(A,B)$ given $S$$$f_{A,B|S}(a,b|s)\propto f_A(a)f_B(s-a)\mathbb{I}_{a+b=s}$$Generating a realisation from this target can be done directly if the shape is simple enough, or by accept-reject, Metropolis-Hastings, slice sampling, or any other standard simulation method.

Answer 3

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.

Answer 4

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