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

def answer(x,y,n):
    n = float(n)
    ans = []
    h = []
    counts1 = []
    counts2 = []
    perms = []

    for i in xrange(int(n)):

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

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

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

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

#         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))

    answer1(x,y,n, c1)

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

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

def answer1(x,y,n,c1):
    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.

  • 3
    $\begingroup$ 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. $\endgroup$
    – psarka
    Dec 5 '14 at 8:57
  • $\begingroup$ 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. $\endgroup$
    – Bill
    Dec 5 '14 at 17:11
  • 1
    $\begingroup$ 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. $\endgroup$
    – 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.

Visibility examples

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.


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.


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:

Distribution for n=40

  • $\begingroup$ 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)? $\endgroup$
    – Vaughn
    Dec 6 '14 at 20:23
  • $\begingroup$ 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.) $\endgroup$
    – whuber
    Dec 7 '14 at 2:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.