# PDF of a sum of dependent variables

This is a direct continuation of my recent question. The thing that I actually want to get is the distribution of $a+d+\sqrt{(a-d)^2+4bc}$, where $a,b,c,d$ are uniform in $[0,1]$. Now, the distribution of $(a-d)^2+4bc$ was successfully computed in the mentioned thread, and let's call it $h(x)$. The distribution of $\sqrt{(a-d)^2+4bc}$ is simply $h(x^2)\cdot 2x$. The last step would be to compute the distribution of the sum of $X=a+d$ and $Y=\sqrt{(a-d)^2+4bc}$ in a way similar to the previous one, but $X$ and $Y$ are not independent, and now I'm stuck and don't even know where to start with.

It may be useful to note that $\sqrt{(a-d)^2+4bc}=\sqrt{(a+d)^2-4(ad-bc)}$ and in the latter the components under the root (i.e., $X^2=(a+d)^2$ and $W=-4(ad-bc)$) are easy to calculate. Then, I'm interested in the distribution of $X+\sqrt{X^2+W}$, knowing the distributions of $X$ and $\sqrt{X^2+W}$.

I don't see any useful change of variables. I thought about using conditional probability, but how can I find $f(\sqrt{X^2+W}\Big|X)$? I might be too much ahead and maybe have to go back a few steps.

Is it even possible to calculate something like this?

The resulting distribution should look like this:

EDIT: The accepted answer gives the solution I was looking for, however, I'm still curious how to derive it analytically. I mean, in my previous question the CDF was given as an integral:

$\int_0^4 F(\delta-y)g(y)dy$

with $F$ and $g$ given by simple functions. Theoretically, that could be integrated using pen and paper. Of course using software is natural. However, I'm still curious how to give a closed-form answer here. wolfies answer rings a bell, but... A convolution of three pdfs of such a (relatively) complicated function?

Find the pdf of: $\quad A+D+\sqrt{(A-D)^2+4BC}, \quad$ where $A,B,C,D$ are iid $Uniform(0,1)$

Let $U = 4 BC$, where $U$ has pdf: $\quad g(u) = \frac{1}{4} \log \left(\frac{4}{u}\right) \quad \text{for } 0<u<4$.

This reduces the problem from 4 to 3 independent random variables. Then, by independence, the joint pdf of $(A,D,U)$ is $f(a,d,u)$:

Let $Z = A+D+\sqrt{(A-D)^2+4BC}$. The cdf of $Z$ is $P(Z<z)$:

where I am using the Prob function from the mathStatica package for Mathematica to automate the nitty-gritties.

The pdf of $Z$ is simply the derivative of the latter wrt $z$, which yields the solution:

All done.

Here is a plot of the exact theoretical pdf of $Z$:

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.

• Neat! Although I don't have mathStatica, I managed to do it directly in Mathematica. This answers my question rather fully, but I'm still curious how to do it without a computer, in a manner similar to my previous question. There, whuber gave the integral explicitly and theoretically that could be calculated using a pen and paper. Of course using software is natural, but how should I proceed in the current case? Apr 1, 2015 at 16:19
• Abramowitz and Stegun? ;) Apr 1, 2015 at 16:35
• Well, you know the integrals that need to be evaluated, so it is just a question of evaluating them. In the days before we had computer algebra systems, when faced with unpleasant and tricky integration tasks outside of the ordinary, one typically headed off to tables of integrals such as Abramowitz and Stegun. Apr 1, 2015 at 17:18

Just after reading wolfies answer I understood I could calculate the final distribution from the very beginning without all the mid-point steps:

M[x_] := M[x] = Evaluate@FullSimplify@ Integrate[ Boole[a + d + Sqrt[(a - d)^2 + 4 b c] <= x], {a, 0, 1}, {b, 0, 1}, {c, 0, 1}, {d, 0, 1}] gives the CDF and

m[x_] := m[x] = Evaluate@FullSimplify@D[M[x], x] gives the PDF that works perfect with my simulation:

This uses directly the approach of an answer to my previous question.

• Yes - that works nicely here. But interestingly, it does NOT seem to work for your original (simpler) problem. That is, Integrate[ Boole[(a-d)^2 + 4 b c < x], {a,0,1}, {b,0,1}, {c,0,1}, {d,0,1}] returns an unevaluated integral. Apr 1, 2015 at 16:32