Let $X_i$, $i=1,2,3$ be three independent exponentially distributed random variables with parameter $\lambda_i$. Define $R_1=\min\{X_1,X_3\}$ and $R_2=\min\{X_2,X_3\}$.
Then, $R_1$ and $R_2$ are obviously not independent. The distribution function is: $$P(R_1<r_1,R_2<r_2)=1-P(R_1>r_1)-P(R_2>r_2)+P(R_1>r_1,R_2>r_2),$$ where $P(R_i>r_i)=e^{-(\lambda_i+\lambda_3)r_i}$, $i=1,2$ and $P(R_1>r_1,R_2>r_2)=e^{-\lambda_1 r_1}e^{-\lambda_2 r_2}e^{-\lambda_3 \max\{r_1,r_2\}}$.
Now, this means that the distribution function $F_{(R_1,R_2)}(r_1,r_2)$ is given by $$F_{(R_1,R_2)}(r_1,r_2)=\begin{cases} 1 -e^{-(\lambda_1+\lambda_3)r_1}-e^{-(\lambda_2+\lambda_3)r_2}+ e^{-\lambda_1 r_1}e^{-(\lambda_2+\lambda_3) r_2}, \mbox{ if } r_1<r_2 \\ 1-e^{-(\lambda_1+\lambda_3)r_1}-e^{-(\lambda_2+\lambda_3)r_2}+e^{-(\lambda_1+\lambda_3) r_1}e^{-\lambda_2 r_2}, \mbox{ if } r_1>r_2 \end{cases}.$$
However, when I differentiate this function to obtain the density, i.e $f_{(R_1,R_2)}(r_1,r_2)=\frac{\partial^2}{\partial r_1 \partial r_2}F_{(R_1,R_2)}(r_1,r_2)$ I do not get a function that integrates to $1$ on $(0,\infty)^2$. I get, $$f_{(R_1,R_2)}(r_1,r_2)=\begin{cases} \lambda_1 (\lambda_2+\lambda_3)e^{-\lambda_1 r_1}e^{-(\lambda_2+\lambda_3) r_2}, \mbox{ if } r_1<r_2 \\ \lambda_2(\lambda_1+\lambda_3)e^{-(\lambda_1+\lambda_3) r_1}e^{-\lambda_2 r_2}, \mbox{ if } r_1>r_2 \end{cases}.$$
Observe that $F$ is indeed a distribution function (continuous, increasing, $F(0,r_2)=F(r_1,0)=0$ and $F(r_1,r_2)\to 1$ as $(r_1,r_2)\to (\infty,\infty)$. So I expect that the derivative should give me a proper density function, but I get $$\int_0^\infty \int_0^\infty f_{(R_1,R_2)}(r_1,r_2) dr_1 dr_1= 1-\frac{\lambda_3}{\lambda_1+\lambda_2+\lambda_3}.$$
Where is the mistake? I must admit I'm a bit puzzled about this... there must be a very silly mistake somewhere. Thanks!
Are you sure there is a density function?
Some of the time $X_3$ will be the smallest of the three $X_i$, and when that happens $R_1=R_2$. That is, $P[R_1=R_2]>0$, and some of the probability mass of the joint distribution of $(R_1,R_2)$ is concentrated on the set $\{(r,r):r>0\}$, a set of Lebesgue measure zero in the plane.
If you have the visualization tools handy, you can generate a thousand or so Monte Carlo sample values of $(R_1,R_2)$ and plot a 2-D histogram to see what the empirical density looks like.