What's the distribution of $(a-d)^2+4bc$, where $a,b,c,d$ are uniform distributions?

Question

I have four independent uniformly distributed variables $a,b,c,d$, each in $[0,1]$. I want to calculate the distribution of $(a-d)^2+4bc$. I computed the distribution of $u_2=4bc$ to be $$f_2(u_2)=-\frac{1}{4}\ln\frac{u_2}{4}$$ (hence $u_2\in(0,4]$), and of $u_1=(a-d)^2$ to be $$f_1(u_1)=\frac{1-\sqrt{u_1}}{\sqrt{u_1}}.$$ Now, the distribution of a sum $u_1+u_2$ is ($u_1,\, u_2$ are also independent) $$f_{u_1+u_2}(x)=\int_{-\infty}^{+\infty}f_1(x-y)f_2(y)dy=-\frac{1}{4}\int_0^4\frac{1-\sqrt{x-y}}{\sqrt{x-y}}\cdot\ln\frac{y}{4}dy,$$ because $y\in(0,4]$. Here, it has to be $x>y$ so the integral is equal to $$f_{u_1+u_2}(x)=-\frac{1}{4}\int_0^{x}\frac{1-\sqrt{x-y}}{\sqrt{x-y}}\cdot\ln\frac{y}{4}dy.$$ Now I insert it to Mathematica and get that $$f_{u_1+u_2}(x)=\frac{1}{4}\left[-x+x\ln\frac{x}{4}-2\sqrt{x}\left(-2+\ln x\right)\right].$$

I made four independent sets $a,b,c,d$ consisting of $10^6$ numbers each and drew a histogram of $(a-d)^2+4bc$:

and drew a plot of $f_{u_1+u_2}(x)$:

Generally, the plot is similar to the histogram, but on the interval $(0,5)$ most of it is negative (the root is at 2.27034). And the integral of the positive part is $\approx 0.77$.

Where's the mistake? Or where am I missing something?

EDIT: I scaled the histogram to show the PDF.

EDIT 2: I think I know where's the problem in my reasoning - in the integration limits. Because $y\in (0,4]$ and $x-y\in(0,1]$, I cannot simply $\int_0^x$. The plot shows the region I have to integrate in:

This means I have $\int_0^x$ for $y\in(0,1]$ (that's why part of my $f$ was correct), $\int_{x-1}^x$ in $y\in(1,4]$, and $\int_{x-1}^4$ in $y\in (4,5]$. Unfortunately, Mathematica fails to compute the latter two integrals (well, it does calculate the second, by there's an imaginary unit in the output that spoils everything...).

EDIT 3: It appears that Mathematica CAN compute the last three integrals with the following code:

(1/4)*Integrate[((1-Sqrt[u1-u2])*Log[4/u2])/Sqrt[u1-u2],{u2,0,u1}, Assumptions ->0 <= u2 <= u1 && u1 > 0]

(1/4)*Integrate[((1-Sqrt[u1-u2])*Log[4/u2])/Sqrt[u1-u2],{u2,u1-1,u1}, Assumptions -> 1 <= u2 <= 3 && u1 > 0]

(1/4)*Integrate[((1-Sqrt[u1-u2])*Log[4/u2])/Sqrt[u1-u2],{u2,u1-1,4}, Assumptions -> 4 <= u2 <= 4 && u1 > 0]

which gives a correct answer :)

Answer 1

Often it helps to use cumulative distribution functions.

First,

$$F(x) = \Pr((a-d)^2 \le x) = \Pr(|a-d| \le \sqrt{x}) = 1 - (1-\sqrt{x})^2 = 2\sqrt{x} - x.$$

Next,

$$G(y) = \Pr(4 b c \le y) = \Pr(b c \le \frac{y}{4}) = \int_0^{y/4} dt + \int_{y/4}^1\frac{y\,dt}{4t} = \frac{y}{4}\left(1 - \log\left(\frac{y}{4}\right)\right).$$

Let $\delta$ range between the smallest ($0$) and largest ($5$) possible values of $(a-d)^2 + 4 b c$. Writing $x=(a-d)^2$ with CDF $F$ and $y=4 b c$ with PDF $g = G^\prime$, we need to compute

$$H(\delta) = \Pr((a-d)^2 + 4 b c \le \delta) = \Pr(x\le \delta-y) = \int_0^4 F(\delta-y)g(y)dy.$$

We can expect this to be nasty--the uniform distribution PDF is discontinuous and thus ought to produce breaks in the definition of $H$--so it is somewhat amazing that Mathematica obtains a closed form (which I will not reproduce here). Differentiating it with respect to $\delta$ gives the desired density. It is defined piecewise within three intervals. In $0 \lt \delta \lt 1$,

$$H^\prime(\delta) = h(\delta) = \frac{1}{8} \left(8 \sqrt{\delta }+\delta (-(2+\log (16)))+2 \left(\delta -2 \sqrt{\delta }\right) \log (\delta )\right).$$

In $1 \lt \delta \lt 4$,

$$h(\delta) = \frac{1}{4} \left(-(\delta +1) \log (\delta -1)+\delta \log (\delta )-4 \sqrt{\delta } \coth ^{-1}\left(\sqrt{\delta }\right)+3+\log (4)\right).$$

And in $4 \lt \delta \lt 5$,

$$\eqalign{ &h(\delta) = \\ &\frac{1}{4}\left(\delta -4 \sqrt{\delta -4}+(\delta +1) \log \left(\frac{4}{\delta -1}\right)+4 \sqrt{\delta } \tanh ^{-1}\left(\frac{\sqrt{(\delta -4) \delta }-\sqrt{\delta }}{\delta -\sqrt{\delta -4}}\right)-1\right). }$$

This figure overlays a plot of $h$ on a histogram of $10^6$ iid realizations of $(a-d)^2 + 4bc$. The two are almost indistinguishable, suggesting the correctness of the formula for $h$.

---

The following is a nearly mindless, brute-force Mathematica solution. It automates practically everything about the calculation. For instance, it will even compute the range of the resulting variable:

ClearAll[ a, b, c, d, ff, gg, hh, g, h, x, y, z, zMin, zMax, assumptions];
assumptions = 0 <= a <= 1 && 0 <= b <= 1 && 0 <= c <= 1 && 0 <= d <= 1; 
zMax = First@Maximize[{(a - d)^2 + 4 b c, assumptions}, {a, b, c, d}];
zMin = First@Minimize[{(a - d)^2 + 4 b c, assumptions}, {a, b, c, d}];

Here is all the integration and differentiation. (Be patient; computing $H$ takes a couple of minutes.)

ff[x_] := Evaluate@FullSimplify@Integrate[Boole[(a - d)^2 <= x], {a, 0, 1}, {d, 0, 1}];
gg[y_] := Evaluate@FullSimplify@Integrate[Boole[4 b c <= y], {b, 0, 1}, {c, 0, 1}];
g[y_]  := Evaluate@FullSimplify@D[gg[y], y];
hh[z_] := Evaluate@FullSimplify@Integrate[ff[-y + z] g[y], {y, 0, 4}, 
          Assumptions -> zMin <= z <= zMax];
h[z_]  :=  Evaluate@FullSimplify@D[hh[z], z];

Finally, a simulation and comparison to the graph of $h$:

x = RandomReal[{0, 1}, {4, 10^6}];
x = (x[[1, All]] - x[[4, All]])^2 + 4 x[[2, All]] x[[3, All]];
Show[Histogram[x, {.1}, "PDF"], 
 Plot[h[z], {z, zMin, zMax}, Exclusions -> {1, 4}], 
 AxesLabel -> {"\[Delta]", "Density"}, BaseStyle -> Medium, 
 Ticks -> {{{0, "0"}, {1, "1"}, {4, "4"}, {5, "5"}}, Automatic}]

Answer 2

Like the OP and whuber, I would use independence to break this up into simpler problems:

Let $X = (a-d)^2$. Then the pdf of $X$, say $f(x)$ is:

Let $Y = 4 b c$. Then the pdf of $Y$, say $g(y)$ is:

The problem reduces to now finding the pdf of $X + Y$. There may be many ways of doing this, but the simplest for me is to use a function called TransformSum from the current developmental version of mathStatica. Unfortunately, this is not available in a public release at the present time, but here is the input:

TransformSum[{f,g}, z]

which returns the pdf of $Z = X + Y$ as the piecewise function:

Here is a plot of the pdf just derived, say $h(z)$:

Quick Monte Carlo check

The following diagram compares an empirical Monte Carlo approximation of the pdf (squiggly blue) to the theoretical pdf derived above (red dashed). Looks fine.