Nested normal-distribution integral

Is there an analytical or approximate solution of the following integral?

$$\int_{-\infty}^{\infty}\int_{y-d}^{y+d}\exp\big(-{(x-\mu_1)^2}/{2\sigma^2}\big) \exp\big(-{(y-\mu_2)^2}/{2\sigma^2}\big) \,\mathrm{d}x\,\mathrm{d}y$$

-
Is $x$ supposed to be in the argument of both of the exponentials? – cardinal May 30 '11 at 1:52
@cardinal: Sorry I made a mistake, in first exponent x is the variable and in the second exponent y – shaikh May 30 '11 at 2:35
No problem. Just wanted to make sure. :) – cardinal May 30 '11 at 2:37
In general, no analytical solution exists. This was proved by Liouville in the early 1800s. – cardinal May 30 '11 at 2:44

You can solve this using the same approach used in this answer.

EDIT 1:

The iterated integral can be expressed simply as $$I = 2\pi \sigma ^2 \bigg[\Phi \bigg(\frac{{d - (\mu _1 - \mu _2 )}}{{\sqrt {2\sigma ^2 } }}\bigg) - \Phi \bigg(\frac{{ - d - (\mu _1 - \mu _2 )}}{{\sqrt {2\sigma ^2 } }}\bigg)\bigg],$$ where $\Phi$ is the distribution function of the standard normal distribution (confirmed numerically). I'll provide the derivation later on.

EDIT 2: Before I give the derivation of this result, here are a few numerical results confirming it. Let $r_1 = r_1(d,\mu_1,\mu_2,\sigma)$ denote the approximation obtained for the iterated integral and $r_2 = r_2(d,\mu_1,\mu_2,\sigma)$ the approximation obtained for the simple expression of $I$ above. The following results were obtained: $$d=1.2, \mu_1=2.8, \mu_2=1.6, \sigma=1.9: r_1 = 7.125001170659572 , r_2 = 7.125000629782468$$ $$d=1.9, \mu_1=2.4, \mu_2=-4.6, \sigma=1.5: r_1 = 0.11438563902950714 , r_2 = 0.11438586127326182$$ $$d = 3.2, \mu_1=-0.4, \mu_2=8.1, \sigma=2.9: r_1 = 5.070682528645151 , r_2 = 5.070686434762869$$ $$d=5.2, \mu_1=-1.4, \mu_2=-3.1,\sigma=0.8: r_1 = 4.017263196283472 , r_2 = 4.017262629040601.$$

EDIT 3:

Proof: First write the iterated integral as $$I = 2\pi \sigma ^2 \int_{ - \infty }^\infty {\bigg[\int_{y - d}^{y + d} {f_X (x)\,dx} \bigg]f_Y (y)\,dy} ,$$ where $f_X$ and $f_Y$ are the probability density functions of the ${\rm N}(\mu_1,\sigma^2)$ and ${\rm N}(\mu_2,\sigma^2)$ distributions, respectively. Hence, $$I = 2\pi \sigma ^2 \int_{ - \infty }^\infty {{\rm P}( - d \le X - y \le d)f_Y (y)\,dy} ,$$ where $X \sim {\rm N}(\mu_1,\sigma_2)$. In turn, by the law of total probability, $$I = 2\pi \sigma ^2 {\rm P}( - d \le X - Y \le d),$$ where $Y \sim {\rm N}(\mu_2,\sigma^2)$ and is independent of $X$. Finally, since $X-Y \sim {\rm N}(\mu_1 - \mu_2,2\sigma^2)$, $$I = 2\pi \sigma ^2 {\rm P}\bigg(\frac{{ - d - (\mu _1 - \mu _2 )}}{{\sqrt {2\sigma ^2 } }} \le Z \le \frac{{d - (\mu _1 - \mu _2 )}}{{\sqrt {2\sigma ^2 } }}\bigg),$$ where $Z$ is a standard normal random variable. The result is thus established.

-
Thanks a lot Mr Shai Covo, can you please explain me how you have obtained the above expression out of iterated integral? The first exponent is a function of x whereas the second exponent is the function of y, did you notice that? – shaikh May 30 '11 at 4:16
shaikh: Yes, I noticed that; I'll give the details soon. – Shai Covo May 30 '11 at 5:13
@Covo: it seems quite accurate. I am desperately waiting for the derivation... :) – shaikh May 30 '11 at 6:36
@Covo: Thanks a lot for your help. I could not get your 3rd step of proof which includes the total law of probability. How can $$I = 2\pi \sigma ^2 \int_{ - \infty }^\infty {{\rm P}( - d \le X - y \le d)f_Y (y)\,dy}$$ be transferred to $$I = 2\pi \sigma ^2 {\rm P}( - d \le X - Y \le d)$$ – shaikh May 30 '11 at 7:36
@shaikh: Later on today (but not soon), I will add some details concerning the law of total probability. – Shai Covo May 31 '11 at 13:39

In response to OP's comment, I elaborate here on the equality $${\rm P}( - d \le X - Y \le d) = \int_{ - \infty }^\infty {{\rm P}( - d \le X - y \le d)f_Y (y)\,dy}.$$ A particular case of the law of total probability allows us to compute the probability of event by conditioning on a suitable random variable, as follows: $${\rm P}(A) = \int_{ - \infty }^\infty {{\rm P}(A|Z = z)f_Z (z)\,dz} ,$$ where $f_Z$ is the probability density function of $Z$. (See also p. 1 here; the law of total probability is a special case of the law of total expectation.) Now, letting $A$ be the event that $-d \leq X-Y \leq d$ and letting $Z=Y$, the above equation yields $${\rm P}( - d \le X - Y \le d) = \int_{ - \infty }^\infty {{\rm P}(- d \le X - Y \le d|Y = y)f_Y (y)\,dy}.$$ Since $X$ and $Y$ are independent, $${\rm P}( - d \le X - Y \le d|Y = y) = {\rm P}( - d \le X - y \le d).$$ The desired equality is thus established.

-
@Covo: I must say great. A lot of thanks for your time and efforts. I got it now. Again thanks for providing all these explanations. – shaikh Jun 1 '11 at 1:57
@shaikh: Thanks, glad to help! – Shai Covo Jun 1 '11 at 6:50
@Covo: How did you find the approximated numerical results for the iterated integral in your first answer Edit2? Which approximation did you use? – shaikh Jun 1 '11 at 9:52
@shaikh: I have programmed the approximations in Java. If we define $\varphi(y)=\varphi(y;d,\mu_1,\sigma)$ by $\varphi (y) = \int_{y - d}^{y + d} {\exp ( - \frac{{(x - \mu _1 )^2 }}{{2\sigma ^2 }})dx}$, then the iterated integral becomes $\int_{ - \infty }^\infty {\varphi (y)\exp ( - \frac{{(y - \mu _2 )^2 }}{{2\sigma ^2 }})dy}$. So, I programmed a function that returns an approximated value for $\varphi(y)$, thus reducing the problem to approximating an ordinary improper integral from $-\infty$ to $\infty$; as the limits of integration I took $\mp 100$ (instead of $\mp \infty$). – Shai Covo Jun 1 '11 at 10:26
public static double getPhi(double y, double d, double mu1, double sigma){ double s=0; int N=500; double h=(2*d)/N; double x=0; for(int i=0; i<N; i++){ x=(y-d)+(i+0.5)*h; s=s+Math.exp(-((x-mu1)*(x-mu1))/(2*sigma*sigma))*h; } return(s); } – Shai Covo Jun 1 '11 at 12:16