Descriptive statistics of two events

Question

I am trying to use the descriptive statistics feature of Mathematica to answer the following question: suppose I have two events, A and B, whose occurrence is described by a normal distribution around time $t_A$ and $t_B$ (with moments $\sigma_A$ and $\sigma_B$). What is the probability that, at a given time $t$, both events have occurred?

I have tried to define distributions and then use them in a construct of the form Probability[A && B, …], but it is clear from the documentation that Probability is not supposed to be used that way.

---

PS: I'm not so much interested in the mathematical solution (I have to read up on that subject, I'm really bad in statistics), rather in the way I can get Mathematica to solve this issue without me doing the maths :)

Answer 1

For reasons @wolfies has amply made clear, it's wise to solve probability problems in Mathematica by getting as close to first principles as possible, rather than relying on black-box solutions (which tend to be much slower as well).

Assuming (as the O.P. implicitly does) that the events are independent, the axioms of probability assert that the chance both events have occurred is the product of their chances of occurring. By definition, the chance that an event has occurred by time $t$ is the value of its cumulative distribution function at $t$. Therefore, after consulting the help page to make sure about the parameterization of Normal distributions, we can reliably construct a correct answer as

f[t_, {tA_, sA_}, {tB_, sB_}] :=
   Evaluate[CDF[NormalDistribution[tA, sA], t] CDF[NormalDistribution[tB, sB], t]];

(I stuck in the Evaluate so we can see what MMA's final formula might be and compare it to other putative solutions:

? f

$f[\text{t$\_$},\{\text{tA$\_$},\text{sA$\_$}\},\{\text{tB$\_$},\text{sB$\_$}\}]\text{:=}\frac{1}{4} \text{erfc}\left(\frac{\text{tA}-t}{\sqrt{2} \text{sA}}\right) \text{erfc}\left(\frac{\text{tB}-t}{\sqrt{2} \text{sB}}\right)$

Notice the lack of square roots over $s_A$ and $s_B$.)

To help understand it, let's plot (a) this function against time $t$ and (b) the two probability density functions for the events, allowing manipulation of the four parameters $t_A, \ldots, s_B$. The principal contribution of this code is to determine reasonable ranges for the plots automatically:

Manipulate[
 range = {t, Min[tA, tB] - 3 Sqrt[sA^2 + sB^2], Max[tA, tB] + 3 Sqrt[sA^2 + sB^2]};
 GraphicsRow[{
   Plot[f[t, {tA, sA}, {tB, sB}], Evaluate@range, PlotStyle -> Thick, 
    AxesLabel -> {"Time", "Probability"}, ImageSize -> i],
   Plot[{PDF[NormalDistribution[tA, sA]][t], 
     PDF[NormalDistribution[tB, sB]][t]}, Evaluate@range, 
    Filling -> Axis, PlotRange -> {Full, Full}, 
    AxesLabel -> {"Time", "Density"}, ImageSize -> i]
   }],
 {{tA, 0}, -5, 0}, {{tB, 2}, 0, 5}, {{sA, 3}, 0, 3}, {{sB, 1/2}, 0, 3}, 
   {{i, 300, "Image size"}, 50, 500}
 ]

The illustration makes sense: there is essentially no probability of both events occurring before time $t=1$, when the event $B$ (red) first begins to have some visible chance of happening. At that point the chance of both events (left) rises rapidly because it's likely event $A$ (blue) has already occurred and the chance of event $B$ is rapidly increasing (as attested by the height of its density curve). By time $t=4$, the chance is good both events have already occurred, with some residual uncertainty about $A$ due to its wide spread. The righthand density plots appear correctly to represent Normal distributions with the given locations ($0$ and $2$) and spreads ($3$ and $1/2$, respectively) while the foregoing reasoning indicates the lefthand plot--our solution--correctly reflects the intended result.

Answer 2

The above answer by chris is unfortunately somewhat confused. I would prefer to make this as a comment, but the comment box does not allow for inputs and outputs, so this is the best I can do. I do not think it fair to apportion blame to the author, for the notation he has used is, in fact, entirely natural. It is also the notation used by the second answer! Rather, the fault lies with the Mathematica syntax for these distributions which is frankly appalling and internally inconsistent in its design. The author above (chris) is just another victim of that design.

In particular, the OP considers two Normal random variables, say X and Y, which Chris enters into Mathematica as:

AA = PDF[NormalDistribution[t1, s1], x]
BB = PDF[NormalDistribution[t2, s2], y]

where, in Mathematica notation, s1 and s2 denote the standard deviation of the respective distributions.

He then considers a bivariate Normal pdf, which he enters, quite naturally, as:

CC = PDF[MultinormalDistribution[{t1, t2}, {{s1, 0}, {0, s2}}], {x,y}]

In the case of independence, this is meant to have the property that CC = AA * BB. But due to a notation kerfuffle in the Mathematica implementation, this is not the case, and Chris' example has gone awry. The following should simplify to 1:

To get what he seeks, one would need to enter the bivariate Normal as:

CCFixed = PDF[MultinormalDistribution[{t1, t2}, {{s1^2, 0}, {0, s2^2}}], {x,y}]

Then:

CCFixed/(AA*BB) // Simplify

will return an entity that simplifies to 1.

The mistake has occurred for essentially 2 reasons: the first is the internally inconsistent syntax between the Normal and Multinormal in mma; the second is the use of 'black boxes' to denote distributional forms which can encourage people to use them without checking the actual output generated. A dangerous combo, as the answers above illustrate.

Answer 3

You can compute this probability using Probability as you hinted

prob = Simplify[ Probability[ x > x1 && x < x2 && y > y1 && y < y2, {x, y} \[Distributed] 
MultinormalDistribution[{t1, t2}, {{σ1^2, 0}, {0, σ2^2}}]], 
 Assumptions -> {σ1 > 0, σ2 > 0}];

Let us give it some numerical values

numtest = {x1 -> 1, x2 -> 1.1, y1 -> 1, 
  y2 -> 1.1, σ1 -> 1, σ2 -> 1, t1 -> 1/2, t2 -> 3/2}

prob /. numtest

(* ==> 0.00123564 *)

Now we can check graphically if it makes sense?

Show[
 Plot[{PDF[NormalDistribution[t1, σ1], x], 
     PDF[NormalDistribution[t2, σ2], x]} /. numtest // 
   Evaluate, {x, -4, 4}],
 Plot[{PDF[NormalDistribution[t1, σ1], x], 
     PDF[NormalDistribution[t2, σ2], x]} /. numtest // 
   Evaluate, {x, 1, 1.1}, Filling -> Axis, PlotRange -> {-0.1, 0.5}]]

qualitatively the 2D volume is 0.1^2*0.35^2=0.00125.

EDIT

Thanks to Wolfies for the typo on sigmas (numerical example was ok by chance).

Answer 4

We can translate your question 'What is the probability that, at a given time $t$, both events have occurred?' to the statement $t\geq t_A \wedge t\geq t_B$, which we can feed to Probability:

Assuming[{Subscript[σ,A]>0,Subscript[σ,B]>0},
  Probability[
    t>=Subscript[t,1]\[And]t>=Subscript[t,2],
    {Subscript[t,1],Subscript[t,2]} \[Distributed]
      MultinormalDistribution[
        {Subscript[t,A],Subscript[t,B]},
        {{Subscript[σ,A]^2,0},{0,Subscript[σ,B]^2}}
      ]
  ]
]

resulting in

$$\frac{1}{4} \left(\text{erf}\left(\frac{t-t_A}{\sqrt{2} \sigma _A}\right)+1\right) \left(\text{erf}\left(\frac{t-t_B}{\sqrt{2} \sigma _B}\right)+1\right)$$

Here we have silently assumed that the events $A$ and $B$ are independent.