# Probability of Two Independent Events, Unknown Distribution

I'm trying to deepen my understand of probability, unions and intersections.

How would you go about finding the probability of x and y ( P(x U y) ??) for a collection of n objects (2 < n < 40) of unique heights, arranged along a line sequentially. Where x is the number of objects visible if you were to look from one direction along the line, and y is the number visible from the opposite direction?

It's a pretty good brain teaser I'm trying to wrap my head around. I wrote some python to help my generate some empirical evidence. I can easily brute force check the number of permutations up to n = 10, but above that permutations are n factorial, which gets large, quickly. So I'm trying to determine the probabilities for X and Y, or even the distinct counts themselves.

import itertools
from collections import Counter
import math

f = math.factorial

n = float(n)
ans = []
h = []
counts1 = []
counts2 = []
perms = []

for i in xrange(int(n)):
h.append(i+1)

for i in itertools.permutations(h,len(h)):
count1 = 1
count2 = 1
perms.append(i)
p1 = i
p2 = i[-1]

for j1,k1 in enumerate(i):
if i[j1] > p1:
p1 = i[j1]
count1 += 1
counts1.append(count1)

rev = [z for z in reversed(i)]

for j2,k2 in enumerate(rev):
if rev[j2] > p2:
p2 = rev[j2]
count2 += 1
counts2.append(count2)

#         print count1, i, count2

#     print zip(counts1, counts2)

c1 = Counter(counts1); print c1, "\n"
cz = zip(counts1, counts2)
czip = Counter(cz); print czip, "\n"
combs = [i for i in itertools.combinations_with_replacement(h,2)]

print "A \tC(A) \tP(A)"
for i in combs:
print i, "\t", czip[i], "\t", czip[i]/float(f(n))

ans = czip[(x,y)]
#     print ans

return ans, float(ans) / f(n)

print "\n##### just calculations #####"

f = math.factorial
fn = float(f(n))
nn = int(fn/n)

print "\n", "n =", n, "\n","n! =", math.factorial(n), "\n", "n!/n =", nn, "\n"

px = c1[x] / fn
print "P(X), x =", x,"\t", c1[x],"\t", px, "\t", px * nn

py = c1[y] / fn
print "P(Y), y =", y,"\t", c1[y],"\t", py, "\t", py * nn

c1xORy = c1[x] + c1[y]
pxORy = c1xORy / fn
print "P(X or Y) \t", c1xORy, "\t", pxORy, "\t", pxORy * nn

c1xANDy = abs(c1[x]-c1[y])
pxANDy =  c1xANDy / fn # px + py - pxORy
print "P(X and Y) \t",c1xANDy, "\t", pxANDy, "\t", pxANDy * nn, "\n"


Insight much appreciated.

• If $X$ is the number of visible objects, then $X$ is a random variable rather than event (and $Y$ as well). There is no meaningful 'union' of random variables, so you may want to replace $P(X\cup Y)$ by something else. Dec 5 '14 at 8:57
• Also, it's going to depend on both how high the eye is above the line and how far apart the objects are, no? I don't read python, so I can't figure out what you are assuming about that.
– Bill
Dec 5 '14 at 17:11
• Note that these random variables (number visible from right and number from the left) are not independent, despite that characterization in the title. They tend to be inversely associated.
– whuber
Dec 5 '14 at 17:37

To make this question definite and answerable, I interpret it as placing $n$ distinct positive numbers $h_1 \lt h_2 \lt \cdots \lt h_n$ independently and randomly on the Real line at locations $x_1, x_2, \ldots, x_n$. I assume the (common) distribution of these random locations is continuous so that there is zero chance of collision. Finally, the number $h_s$ is "visible" from the left provided all the numbers at locations less than $x_s$ are less than $h_s$, with a similar definition of visibility from the right. The figure shows four random configurations of heights $(1, 2, \ldots, 15)$. Those visible from the right are shown in red and those visible from the left in blue. One is always visible from the left and the right; it is shown in black.

### Analysis

Visibility is determined by the sequence in which the heights are encountered from left to right. All sequences are equally likely (regardless of the actual distribution used to place the heights). The chance there are $t$ visible from the left depends on the location of the largest height, so suppose there are $s-1$ heights to its left (that is, the largest height is located at position $s$ counted from the left). Call this chance $e(t, s)$. We count the maximum height as visible and recursively evaluate $e(t,s)$ in terms of the remaining $s-1$ heights to the left of the maximum,

$$e(t,s) = \sum_{i=1}^{s-1} \Pr(t-1\text{ are visible left of } i)\Pr(\text{max height is at }i) \\ = \frac{1}{s-1}\sum_{i=1}^{s-1} e(t-1,i),\, s \gt 1.$$

For $s=1$ it is obvious that $e(1,1)=1$ and $e(t,1)=0$ for $t\ne 1$.

Visibility from both sides can be computed by breaking the configurations into the number visible from the left plus the number visible from the right: they must sum to $t+1$ (because the maximum height is uniquely visible from both sides). This has to be averaged over all the equally likely locations of the maximum. Letting the chance of $t$ heights being visible out of $n$ heights be $p(t,n)$, we conclude

$$p(t,n) = \frac{1}{n}\sum_{s=1}^n \sum_{i=1}^{t+1}e(i,s) e(t+1-i, n+1-s).$$

(This is an example of computing the distribution of a sum of random variables based on their joint distribution.)

The overall computational effort is $O(s^2)$ to find the distribution of $e$ for any given $s$ and therefore $O(n^3)$ to find the distribution of $p$ for a given $n$. This is feasible, taking a few seconds for $n$ up to several hundred or so.

### Implementation

The recursion can be implemented directly in many languages. For those that do not support memoization, it's a good idea to cache the values of $e(t,s)$ as they are computed (effectively making it into a dynamic program). Here is an R example.

e <- function(t, s) {
if (is.na(e.cache[t+1, s+1])) {
e.cache[t+1, s+1] <<- mean(sapply(1:(s-1), function(i) e(t-1,i)))
}
return (e.cache[t+1, s+1])
}

p <- function(t, n) {
mean(sapply(1:n, function(s) {
sum(sapply(1:(t+1), function(i) {
e(i,s) * e(t+1-i, 1+n-s)
}))
}))
}


To use it, create a cache for $e$ and compute $p$:

e.cache <- matrix(NA, n+2, n+2)
e.cache[1:2, ] <- e.cache[, 1:2] <- 0; e.cache[2, 2] <- 1
p.n <- sapply(1:n, function(t) p(t,n))


Here is a partial plot of the distribution stored in p.n: • This is very helpful, thank you @whuber. What are we calculating in p.n? Is it the number visible, given t_x AND t_y (the intersection of t_x and t_y)? Dec 6 '14 at 20:23
• p.n is an array where p.n[t] is the probability that exactly $t$ of the heights will be visible from the ends of the line when there are $n$ heights. ("Intersection" makes no sense here: we are describing a random variable, not a particular event.)
– whuber
Dec 7 '14 at 2:21