Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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

share|cite|improve this question
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
up vote 1 down vote accepted

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


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


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.

share|cite|improve this answer
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.

share|cite|improve this answer
@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

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.