Don't understand the results I'm getting from Integrate

Question

I'm playing around with the standard deviation of the samples taken from the normal distribution. When calculating the integral

$$\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }f(x) f(y) \sqrt{x^2-\frac{1}{2} (x+y)^2+y^2}dydx$$

Where $$f(x)=\text{PDF}[\text{NormalDistribution}[],x]=\frac{e^{-\frac{x^2}{2}}}{\sqrt{2 \pi }}$$

Mathematica returns 0:

f = PDF[NormalDistribution[]];
Integrate[Sqrt[(x^2 + y^2 - (x + y)^2/2)] f[x] f[y], {x, -Infinity, Infinity}, {y, -Infinity, Infinity}]
(*Outputs 0*)

I know that it can't be 0 because for all $x$, $f(x)>0$, and for all $(x,y)\in\mathbb{R}^2$, $\sqrt{x^2-\frac{1}{2} (x+y)^2+y^2} \geq 0$.

If we replace Integrate with NIntegrate, the output is 0.797884, which seems reasonable, but why does Integrate break?

Answer 1

Unfortunately, the problematic part of your integrand is hidden too well to make Mathematica see it. The sqrt-term that resembles a distance from the diagonal is the place where it goes wrong. If you look at it closely and assuming that x and y are real-valued, we get the following simplification

FullSimplify[Sqrt[x^2 + y^2 - 1/2 (x + y)^2], Element[{x, y}, Reals]]
(* Abs[x - y]/Sqrt[2] *)

This term needs to be treated with care. One way, as pointed out by Jim is to use the PrincipalValue option. However, just knowing the fact where the problem lies let's us solve this on our own.

Firstly, we can simplify your integrand and throw it into Mathematica's face that it needs to be careful (my version is 11.2 btw):

f = PDF[NormalDistribution[]];
g = Sqrt[(x^2 + y^2 - (x + y)^2/2)] f[x] f[y];

Integrate[FullSimplify[g, Element[{x, y}, Reals]],
  {x, -Infinity, Infinity}, {y, -Infinity, Infinity}]
(* Sqrt[2/π] *)

Secondly and highly underappreciated is to see the symmetry and integrate only in the valid half of the plane. This result needs than to be multiplied by two:

2*Integrate[g, {x, -Infinity, Infinity}, {y, x, Infinity}] // Simplify
(* Sqrt[2/π] *)

As Geoge pointed out, there is a way of doing the above but calculating both sides:

Integrate[g, {x, -Infinity, Infinity}, {y, -Infinity, x, Infinity}]