Arithmetic with Uncertainty

If you know the temperature in Fahrenheit, then it's easy enough to convert it to Celsius:

def f_to_c(f):
    """Convert Fahrenheit to Celsius."""
    return (f - 32) * 5 / 9

So if it's 68 degrees Fahrenheit, that's the same as f_to_c(68) = 20 degrees Celsius.

But what if we don't know the exact temperature? Maybe we have two thermometers and they disagree by a few degrees: one is showing 66 and the other 71. So we have a rough idea of what the temperature is like in Fahrenheit (somewhere between 66 and 71 degrees), and we'd like to convert this to an equivalently rough idea of what the temperature is like in Celsius.

In a sense, we'd like to simultaneously consider every possible value in the range [66, 71], calculate the equivalent Celsius temperature, and form a range from the resulting values. But because f_to_c is monotonic, we can simply compute the values at the endpoints and use those:

So knowing that the temperature is somewhere between 66 and 71 degrees Fahrenheit is the same as knowing that it's between f_to_c(66) ≈ 19 and f_to_c(71) ≈ 22 degrees Celsius.

One way to think about this computation is as a propagation of our uncertainty. Instead of calculating with numbers, like 68, we're operating on ranges of values, like [66, 71]. This is a crude representation of uncertainty to be sure, but — as we'll see — it's enough to begin to illustrate some neat ideas.

Arithmetic with Uncertainty

Our plan is to write an Ind class that represents uncertain or indeterminate values as ranges of numbers:

class Ind:
    """A representation of an uncertain or indeterminate value.

    Ind(a, b) represents our _knowledge_ that some value lies between a and b."""

    def __init__(self, lo, hi):
        self.lo = lo
        self.hi = hi

    def __eq__(self, other):
        a, b = self.lo, self.hi
        c, d = other.lo, other.hi
        return a == c and b == d

    def __str__(self):
        return f"[{self.lo}, {self.hi}]"

    def __repr__(self):
        return f"Ind({self.lo}, {self.hi})"

so that functions like f_to_c "just work" when applied to Ind objects:

f_temp = Ind(66, 71)
print(f_to_c(f_temp))
# => [18.88888888888889, 21.666666666666668]

In order to make this happen, operators like + and * need to "propagate uncertainty" like we've illustrated above. Now an object like Ind(a, b) represents an uncertain value that's somewhere between a and b: it might be as low as a or as high as b, and that's all we know. So given two Inds — [a, b] and [c, d] — their sum might be as small as a + c or as high as b + d:

class Ind:
    # ...

    def __add__(self, other):
        if isinstance(other, Ind):
            a, b = self.lo, self.hi
            c, d = other.lo, other.hi
            return Ind(a + c, b + d)
        else:
            # ...

What about the else case, where we're adding a number to an Ind? As an example, if you're heading to a class that you know will take exactly 50 minutes, and then walking with a friend for somewhere between 30 and 45 minutes, then between class and walking, you'll spend [50 + 30, 50 + 45] = [80, 95]minutes. So adding a number to anInd just "shifts" it by the number's value. One way to implement this is:

def __add__(self, other):
    if isinstance(other, Ind):
        # ...
    else:
        return Ind(self.lo + other, self.hi + other)

But we'll take a slightly different approach here that nonetheless produces the same result: we'll think of numbers — like the 50 minutes in the example above — as certain knowledge, and represent them as Inds for which lo == hi.

class Ind:
    # ...

    def certain(x):
        """Construct a representation of certainty."""
        return Ind(x, x)

Adding an Ind to a number is then a simple matter of promoting the number to a certain Ind and then adding them together:

def __add__(self, other):
    if isinstance(other, Ind):
        # ...
    else:
        return self + Ind.certain(other)

While we're at it, let's implement __radd__ so that we can add numbers to Inds as well:

def __radd__(self, other):
    return Ind.certain(other) + self

At this point, Ind is functional enough to calculate the "class, then walking" example:

50 + Ind(30, 45)
# => [80, 95]

Let's move on to multiplication. Multiplying two intervals is slightly more involved than addition, because the minimum and maximum values of the product of two intervals isn't always the product of the minimum and maximum values. Here are some examples showing several of the many different possible cases:

# [a, b] * [ c,  d] == [ e,  f]
  [1, 2] * [ 3,  5] == [ 3, 10]  # e = a*b  f = b*d
  [2, 4] * [-2, -1] == [-8, -2]  # e = b*c  f = a*d
  [1, 3] * [-1,  1] == [-3,  3]  # e = b*c  f = b*d
# Many more...

So which components form the min and max depends on the signs of the values and their positions relative to 0. The good news is that we don't actually need to decode the pattern here: instead we can just calculate the four possible products and use the smallest for the min and the largest for the max:

class Ind:
    # ...

    def __mul__(self, other):
        if isinstance(other, Ind):
            a, b = self.lo, self.hi
            c, d = other.lo, other.hi
            # Just consider all possible products.
            ac, ad, bc, bd = a * c, a * d, b * c, b * d
            return Ind(min(ac, ad, bc, bd), max(ac, ad, bc, bd))
        else:
            return self * Ind.certain(other)

    def __rmul__(self, other):
        return Ind.certain(other) * self

Let's quickly dispatch some more operations:

class Ind:
    # ...

    def __sub__(self, other):
        return self + -other

    def __neg__(self):
        return Ind(-self.hi, -self.lo)

    def __rsub__(self, other):
        return Ind.certain(other) - self

    def __truediv__(self, other):
        if isinstance(other, Ind):
            return self * other.reciprocal()
        else:
            return self / Ind.certain(other)

    def reciprocal(self):
        a, b = self.lo, self.hi
        a_inv, b_inv = 1 / a, 1 / b
        return Ind(min(a_inv, b_inv), max(a_inv, b_inv))

    def __rtruediv__(self, other):
        return Ind.certain(other) / self

These are fairly straightforward: we "rewrite" subtraction using addition and negation, and do the same with division in terms of multiplication and reciprocation.

The final operation we'll support is raising to an integer power. This one is slightly tricky: if the power is odd, then this is a monotonic operation, but if it's even, we need to be careful about intervals that straddle 0. As an example, consider the expression [-1, 1] ** 2, which we'll think of as simultaneously squaring every value between -1 and 1. The smallest resulting value is 0, not -1 ** 2:

class Ind:
    def __pow__(self, n):
        a, b = self.lo, self.hi
        an, bn = a**n, b**n

        if n % 2 == 0:
            if a < 0 and b > 0:
                return Ind(0, max(an, bn))
            return Ind(min(an, bn), max(an, bn))
        else:
            return Ind(an, bn)

Here's the complete implementation.

Let's try it out:

f_to_c(Ind(66, 71))
# => [18.88888888888889, 21.666666666666668]

Here are a few more examples: suppose we know that a cylinder's radius is between 10 and 10.5 cm, and that its height is between 12.5 and 13.5 cm. What do we know about its volume?

def cylinder_vol(r, h):
    """Compute the volume of a cylinder with radius r and height h."""
    return 3.14 * r**2 * h


r = Ind(10, 10.5)
h = Ind(12, 13.5)
print(cylinder_vol(r, h))
# => [3768.0, 4673.4975]

Or another: two objects in orbit are between 100 and 102 meters apart, but are drifting towards each other at a rate between 1 and 1.2 meters per second. How long before they collide?

def collision_time(dist, speed):
    """Calculate the time before two objects collide, assuming they're dist m
    apart, and approaching each other at speed m/s."""
    return dist / speed


dist = Ind(100, 102)
speed = Ind(1, 1.2)
print(collision_time(dist, speed))
# => [83.33333333333334, 102.0]

Now each of these problems would have been easy enough to solve manually, i.e. without using our Ind class. In fact, it's clear that all Ind is doing is mildly clever accounting. But the point is that Ind allows us to think at a higher, more abstract level. This is similar to how operations on, say, 3D vectors are "just" operations on their components, but the calculation with vectors is much simpler.

Dependence

Unfortunately there's a subtle problem with our Ind class, related to an issue not present in any of our familiar number systems:

x = Ind(-1, 1)

print(x**2)
# => [0, 1]

print(x * x)
# => [-1, 1]

Squaring a number isn't the same as multiplying it by itself!

Our Ind class isn't smart enough to know that both terms in the expression x * x are the same thing. There's a dependence here that we're not taking into account: the way we compute products assumes that both factors can vary independently of one another. In this case that amounts to assuming that the lefthand factor (x) can be -1 while — at the same time — the righthand factor (also x) can be 1. This isn't possible! But Ind isn't clever enough to detect this.

This notion of independence (and dependence) is a big theme in probability. When calculating with numbers, we don't need to concern ourselves with questions of "identity". But that changes when we consider ranges of values, and will continue to be the case as we consider more exotic, richer models of uncertainty.

Notice how x**2 is a "tighter" interval than x * x. Tighter intervals represent less uncertainty. This makes sense given that __pow__ has more information available to it than __mul__ (namely, that the thing being multiplied by itself n times is the same thing).

Unfortunately, there's no simple fix. One thing we might try is giving each Ind a unique ID, and then check these IDs when operating on them: if the IDs match, we treat them as the same, otherwise we assume that they're independent. This works for very simple calculations like x * x, but it doesn't scale to larger expressions, such as x * y * x + 3 * x. One problem is that compound expressions also have identities (e.g. in x * y + x * y, the two terms x * y are the same); another is that expressions might have "partial dependencies" on each other, like x * y and x.

Losing Information

While we're at it, here's another shortcoming: suppose you're waiting in line at the DMV. There are between 16 and 18 people in front of you (you're not sure if a few people are just there for moral support), and it appears to be taking between 2 and 15 minutes for each person to be served. How long are you going to be waiting?

def wait_time(line_len, wait_per_person):
    """Calculate how long you'll be waiting in line."""
    return line_len * wait_per_person


line_len = Ind(16, 18)
wait = Ind(2, 15)
print(wait_time(line_len, wait))
# => [32, 270]

270 minutes! You might be here a while...

But you probably won't be. The problem is that we're only going to be waiting 270 minutes if our most pessimistic estimations are correct: the line is as long as possible, and everyone takes 15 minutes. We probably don't think everyone will take 15 minutes, but rather that most people will take around 4-5 minutes, with a few people being served extra quickly, and a few very slowly. So our beliefs often have more of a structure to them, a kind of "shape" than we're unable to represent with Ind.

Ind is just too crude for all but the simplest problems. Nonetheless it illustrates several important ideas in probability:

• The value of defining operations on compound or complex objects, and how these can be used to "propagate uncertainty" through transformations.
• The problems posed by dependencies between the quantities we're operating on: we need to consider how the quantities vary with respect to one another.
• The need for a richer representation of uncertainty, which will turn out to involve distributions.

Addendum

Our Ind class implements interval arithmetic, a well-known technique in numerical analysis. The dependency problem described above really is a tough nut to crack.