# Mathematical software for calculations with random variables

I often need to do calculations such as the following:

• $A, B, C, D$ are independent random variables distributed uniformly in [0,1].
• $X := A+B$ and $Y:=B+C+D$.
• What is the probability that $X<Y$?

Is there a free mathematical software that can do such calculations easily?

My favorite mathematical software is SageMath, but as far as I know, its abilities in probability theory are not sufficient for such calculations.

• I guess any software around has a 'random number' generator as $\texttt{Math.random()}$, $\texttt{rand()}$, $\texttt{RandomReal[]}$, etc$\ldots$ – Felix Marin Oct 9 '16 at 21:08
• @FelixMarin I do not need a random number generator... I need to calculate probabilities analytically. – Erel Segal-Halevi Oct 10 '16 at 5:24
• Sorry. I misunderstood your post. Thanks. – Felix Marin Oct 10 '16 at 6:08

If you would like to do these calculations numerically, I think R is very good (everything free).

Here is a link: https://www.r-project.org/

Here is an editor you may like: https://www.rstudio.com/

A solution in R could look like the following

n <- 1000

A <- runif(n, min = 0, max = 1)
B <- runif(n, min = 0, max = 1)
C <- runif(n, min = 0, max = 1)
D <- runif(n, min = 0, max = 1)

X <- A+B
Y <- B+C+D

T <- (X < Y)
sum(T) / length(T) # I get 0.8295

5/6 # True value is 0.8333...


The true value is given by LutzL in the other post.

• This implementation looks clean and clear. – Erel Segal-Halevi Oct 29 '16 at 16:53

$\newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$

You can do it with a simple $\texttt{javascript}$ code which runs in a terminal with $\texttt{node}$ as

$> \texttt{node random0.js}$

// random0.js Felix Marin
"use strict";
var a = null;
var c = null;
var d = null;
var i = 0;
var ITER = 1000000000; // Total number of iterations.
var total = 0;

while (i < ITER) {
a = Math.random();
c = Math.random();
d = Math.random();
if (a < (c + d)) ++total;
++i;
}

console.log("Result = " + total/ITER);

/*
https://math.stackexchange.com/questions/1959638/mathematical-software-for-calculations-with-random-variables
*/

Result = 0.833326178


The result can changes slightly between 'runs' of the script but they remain close to $\bbox[5px,#eee]{\texttt{0.8333}} \approx 5/6$.

In general, you can follow the $\texttt{LutzL}$ 'graphical answer or perform the following evaluation:

\begin{align} &\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\int_{0}^{1} \bracks{a + b < b + c + d}\,\dd d\,\dd c\,\dd b\,\dd a = \int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\bracks{d > a - c}\,\dd d\,\dd c\,\dd a \\[5mm] = &\ \int_{0}^{1}\int_{0}^{1}\braces{% \bracks{a - c < 0}\int_{0}^{1}\,\dd d + \bracks{0 < a - c < 1}\int_{a - c}^{1}\,\dd d}\,\dd c\,\dd a \\[5mm] = &\ \int_{0}^{1}\int_{0}^{1}\bracks{c > a}\,\dd c\,\dd a + \int_{0}^{1}\int_{0}^{1}\bracks{a - 1 < c < a}\pars{1 - a + c}\,\dd c\,\dd a \\[5mm] = &\ \int_{0}^{1}\int_{a}^{1}\,\dd c\,\dd a + \int_{0}^{1}\int_{0}^{a}\pars{1 - a + c}\,\dd c\,\dd a = \int_{0}^{1}\pars{1 - a}\,\dd a + \int_{0}^{1}\pars{a - {1 \over 2}\,a^{2}}\,\dd a \\[5mm] = &\ {1 \over 2} + \pars{{1 \over 2} - {1 \over 6}} = \bbox[5px,border:1px groove navy]{5 \over 6} \end{align}

Set $A'=1-A$ which is equally a uniform random variable on $[0,1]$, then $X<Y$ is equivalent to $$1<A'+C+D.$$

The probability is the volume of the cube minus one edge pyramid. Which comes up to $$1-\frac16=\frac56.$$

• This solves the specific example problem that I gave, but I was looking for a general tool that can automatically solve more complex problems. – Erel Segal-Halevi Oct 29 '16 at 16:52