Probability density function of $\frac{X}{X+Y}$ ($X$ and $Y$ are independently uniformly distributed over $[0,1]$) If $X$ and $Y$ are independent uniform random variables over $[0,1]$, what is the density function of the variable $Z=\frac{X}{X+Y}$? I got the following result:
$$f_Z(z) = 
\begin{cases}
\frac{z^2}{2(1-z)^2}, & \text{if } 0 \le z \le \frac{1}{2} \\
\frac{1}{2}, & \text{if } \frac{1}{2} < z \le 1 \\ 
0, & \text{otherwise}.
\end{cases}$$
My method is as follows: I first write $Z=\frac{1}{1+Y/X}$, and then I calculate the probability function of $\frac{Y}{X}$, and this can be found here. And then I calculate $Z$'s domain and density function according to the domain and density function of $\frac{Y}{X}$ (i.e., let $W=\frac{Y}{X}$, and you can get $W=\frac{1}{Z}-1$, and substitute this into $W$'s domain and density function, you can get $Z$'s). For your convenience, the density function of $W$ is as follows:
$$f_W(w) = 
\begin{cases}
\frac{1}{2}, & \text{if } 0 \le w \le 1 \\
\frac{1}{2w^2}, & \text{if } w > 1 \\ 
0, & \text{otherwise}.
\end{cases}$$
However, when I use Python to generate 100000 samples of $Z$ to test the result, I get the following graph, which contradicts with the density function of $f_Z(z)$ (the $z>0.5$ part). 
The code is as follows:
import numpy as np
from matplotlib import pyplot as plt

a=np.random.rand(100000,2)
s=a.sum(1)
q=a[:,0]/s
plt.hist(q,100)
plt.show()

Where am I wrong? Did I miscalculate the density function of $Z$?
 A: I can't follow your deduction, but the result is wrong.
A simple way is to write the distribution:
$$F_Z(z)=P(Z\le z) = P(X (1-z) \le Y z)$$
Do a graph and realize that that event corresponds to a triangle inside a unit square, compute its area (you have two case: for $z<1/2$  and $z>1/2$) and you should  get
$$P(Z\le z)=\begin{cases}
1- \frac12 \frac{1-z}{z} & z \ge  \frac12 \\
\frac12 \frac{z}{1-z} & z <  \frac12
\end{cases}$$
The derivative of this gives $f_Z(z)$, which corresponds to the graph from your simulation.
Edit: Regarding your derivation. Letting $U=Y/X$, we have $$f_U(u) = 
\begin{cases}
1/2, & \text{if } 0 \le u \le 1 \\
1/(2u^2), & \text{if } u > 1  
\end{cases}$$
Then, for $Z=g(U)=1/(1+U)$ with $|g'(u)|=\frac{1}{(1+u)^2}$ :
$$f_Z(z)= \frac{f_U(u)}{|g'(u)|}=\begin{cases}
(1+u)^2/2, & \text{if } 0 \le u \le 1 \\
(1+u)^2/(2u^2), & \text{if } u > 1  
\end{cases}$$
Replacing $u=(1-z)/z$ you get the density.
A: Your function $p$ does not have integral one, hence it can't be the correct density.
If I'm not mistaken, one should find
$$p(z)=\frac{1}{2\max(z,(1-z))^2}
$$
for any $z\in [0,1]$ (thanks for the corrections in the comments). You should indeed show us your work so we could point out the mistake.
Derivation of the formula
Let $f \colon \mathbb{R} \to \mathbb{R}$ be a bounded measurable function. Then we have (via the substitution $y=x/t-x$, then use of Fubini's theorem)
$$
\begin{align}
\mathbb{E}(f(Z)) &= \int_0^1 \,\mathrm{d}x \int_0^1 f\left(\frac{x}{x+y}\right) \,\mathrm{d}y \\
&= \int_0^1 \,\mathrm{d}x \int_{x/(x+1)}^1 f(t) \left(\frac{x}{t^2} \,\mathrm{d}t\right) \\
&= \int_0^1 x\,\mathrm{d}x \int_{0}^1 \frac{f(t)}{t^2} \mathbb{1}_{\{t>x/(x+1)\}} \,\mathrm{d}t \\
&= \int_0^1 \frac{f(t)}{t^2}\,\mathrm{d}t \int_{0}^1 x \mathbb{1}_{\{x<t/(1-t)\}} \,\mathrm{d}x \\
&= \int_0^1 \frac{f(t)}{t^2}\,\mathrm{d}t \int_{0}^{\min(t/(1-t),1)} x \,\mathrm{d}x \\
&= \int_0^1 \frac{f(t)}{2t^2} \min\left(\frac{t}{1-t},1\right)^2 \,\mathrm{d}t \\
\end{align}
$$
Hence the density is
$$
p(t)=\frac{\min\left(\dfrac{t}{1-t},1\right)^2}{2t^2}=\frac{1}{2\max(t,(1-t))^2}
$$
A: The $z^2$ term should be in the denominator of the other case.
$f_Z(z) = \begin{cases}
\frac{1}{2 (z - 1)^2}, & \text{if } 0 \le z \le \frac{1}{2} \\
\frac{1}{2 z^2}, & \text{if } \frac{1}{2} < z \le 1 \\ 
0, & \text{otherwise}.
\end{cases}$
