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

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

• I like that you've tried checking the reasonableness of your answer by simulation. Your problem is that you know you've made an error, but can't see quite where. Have you considered that you can check each stage of your method, to troubleshoot where the error lies? For example, does the error lie in your $f_1(u_1)$? Well, you can check your calculated PDF against simulated results just like you did for your final answer. Ditto for $f_2$. If $f_1$ and $f_2$ are both correct, then you made the error when combining them. Such step-by-step checking lets you pin-point where you've gone wrong! – Silverfish Mar 30 '15 at 23:50
• I threw away my first attempt and recalculated it from a scratch. I believe $f_1$ and $f_2$ are correct, although I had to manually multiply my initial $f_1$ by 2 to have it normalized to unity. But that just changes the height and does not explain why I have negative $f$. – corey979 Mar 30 '15 at 23:54
• When generating such histograms to compare to calculated algebraic quantities, scale the histogram to a be valid density (and superimpose them if you can). Do a similar check for your f1 and f2 to make sure you have those right; if they're right (I didn't see any good reason to suspect them yet, but its best to double check), then the problem must be later. – Glen_b Mar 30 '15 at 23:55

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}]

• (+1), especially for reminding people that, instead say of density convolutions, "Often it helps to use cumulative distribution functions" -especially when they have such a simple form as here. And you were damn quick, also. – Alecos Papadopoulos Mar 31 '15 at 15:51
• That looks like a neat solution that I'd love to accept - right after I understand it. I'm more a calculus man than a probabilist; at this moment I have three questions: i) how did you use the CDF to get $F(x)$ and $G(y)$, ii) why there's $F$ and $g$ under the integral for $H$, and iii) how do you from its form that the solution result will be piecewise? – corey979 Mar 31 '15 at 15:54
• (1) $F$ and $G$ are the CDFs. They are computed from the definition of a CDF, as indicated by the first equalities following their first appearances. The details should be apparent in the code I have inserted. (2) This is the convolution formula for a sum (more fully explained in a similar calculation at stats.stackexchange.com/a/144237). (3) I inserted a link to another thread about properties of uniform distributions. – whuber Mar 31 '15 at 16:34

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.