Garrison Keillor’s fictional Lake Wobegon is a place “where all the children are above average.” Donald Knuth alluded to this in his exercise regarding “Lake Wobegon Dice,” a set of dice where the roll of each die is (probably) above average.

Let A be a six-sided die with a 5 on one side and 3’s on the other sides.

Let B and C be six-sided dice with 1’s on two sides and 4’s on four sides.

Then the probabilities

Pr( A > (A + B + C)/3 )
Pr( B > (A + B + C)/3 )
Pr( C > (A + B + C)/3 )

are all greater than 1/2.

To see this, list out all the possible rolls where each die is above average. The triples show the value of A, B, and C in that order. We have A above average if the rolls are

(3, 1, 1)
(3, 1, 4)
(3, 4, 1)
(5, *, *)

These outcomes have probability 5/54, 10/54, 10/54, and 1/6. The total is 34/54 = 17/27 > 1/2.

We have B above average when the rolls are

(3, 4, *)
(5, 4, 1)

which have probability 5/9 and 1/27 for a total of 16/27 > 1/2. And C is the same as B.

It seems paradoxical for three different dice to each have a probability better than 1/2 of being above average. The resolution is that the three events

A > (A + B + C)/3
B > (A + B + C)/3
C > (A + B + C)/3

are not exclusive. For example, in the roll (3, 4, 1) both A and B are above the average of 8/3.

I ran across this in The Art of Computer Programming, Volume 4, Pre-fascicle 5A, Exercise 3. Fascicle 5 has since been published, but I don’t have a copy. I don’t know whether this exercise made it into the published version. If it did, the page number may have changed.

More dice posts

• Illustrating Euler characteristic with dice
• More sides or more dice?
• Generating Gaussian samples from dice

A theorem by Paley and Wiener says that a Fourier series with random coefficients produces Brownian motion on [0, 2π]. Specifically,

produces Brownian motion on [0, 2π]. Here the Zs are standard normal (Gaussian) random variables.

Here is a plot of 10 instances of this process.

Here’s the Python code that produced the plots.

    import matplotlib.pyplot as plt
    import numpy as np
    
    def W(t, N):
        z = np.random.normal(0, 1, N)
        s = z[0]*t/(2*np.pi)**0.5
        s += sum(np.sin(0.5*n*t)*z[n]/n for n in range(1, N))*2/np.pi**0.5
        return s
    
    N = 1000
    t = np.linspace(0, 2*np.pi, 100)
    
    for _ in range(10):
        plt.plot(t, W(t, N))
    
    plt.xlabel("$t$")
    plt.ylabel("$W(t)$")
    plt.savefig("random_fourier.png")

Note that you must call the function W(t, N) with a vector t of all the time points you want to see. Each call creates a random Fourier series and samples it at the points given in t. If you were to call W with one time point in each call, each call would be sampling a different Fourier series.

The plots look like Brownian motion. Let’s demonstrate that the axioms for Brownian motion are satisfied. In this post I give three axioms as

1. W(0) = 0.
2. For s > t ≥ 0. W(s) – W(t) has distribution N(0, s – t).
3. For v ≥ u ≥ s ≥ t ≥ 0. W(s) – W(t) and W(v) – W(u) are independent.

The first axiom is obviously satisfied.

To demonstrate the second axiom, let t = 0.3 and s = 0.5, just to pick two arbitrary points. We expect that if we sample W(s) – W(t) a large number of times, the mean of the samples will be near 0 and the sample variance will be around 0.2. Also, we expect the samples should have a normal distribution, and so a quantile-quantile plot would be nearly a straight 45° line.

To demonstrate the third axiom, let u = 1 and v = √7, two arbitrary points in [0, 2π] larger than the first two points we picked. The correlation between our samples from W(s) – W(t) and our samples from W(v) – W(u) should be small.

First we generate our samples.

    N = 1000
    diff1 = np.zeros(N)
    diff2 = np.zeros(N)
    x = [0.3, 0.5, 1, 7**0.5]
    for n in range(N):
        w = W(x)
        diff1[n] = w[1] - w[0]
        diff2[n] = w[3] - w[2]

Now we test axiom 2.

    from statsmodels.api import qqplot
    from scipy.stats import norm

    print(f"diff1 mean = {diff1.mean()}, var = {diff1.var()}.")
    qqplot(diff1, norm(0, 0.2**0.5), line='45')
    plt.savefig("qqplot.png")

When I ran this the mean was 0.0094 and the variance was 0.2017.

Here’s the  Q-Q plot:

And we finish by calculating the correlation between diff1 and diff2. This was 0.0226.

Related posts

• Brownian bridge (interpolating Brownian motion)
• The fractal nature of Brownian motion
• Brown noise and other colors

Yesterday I wrote about how to interpolate a Brownian path. If you’ve created a discrete Brownian path with a certain step size, you can go back and fill in at smaller steps as if you’d generated the values from the beginning.

This means that instead of generating samples in order as a random walk, you could generate the samples recursively. You could generate the first and last points, then fill in a point in the middle, then fill in the points in the middle of each subinterval etc. This is very much like the process of creating a Cantor set or a Koch curve.

In Brownian motion, the variance of the difference in the values at two points in time is proportional to the distance in time. To step forward from your current position by a time t, add a sample from a normal random variable with mean 0 and variance t.

Suppose we want to generate a discrete Brownian path over [0, 128]. We could start by setting W(0) = 0, then setting W(1) to a N(0, 1) sample, then setting W(2) to be W(1) plus a N(0, 1) sample etc.

The following code carries this process out in Python.

    import numpy as np

    N = 128
    W = np.zeros(N+1)
    for i in range(N):
        W[i+1] = W[i] + np.random.normal()

We could make the code more precise, more efficient, but less clear [1], by writing

    W = np.random.normal(0, 1, N+1).cumsum()

In either case, here’s a plot of W.

But if we want to follow the fractal plan mentioned above, we’d set W(0) = 0 and then generate W(128) by adding a N(0, 128) sample. (I’m using the second argument in the normal distribution to be variance. Conventions vary, and in fact NumPy’s random.gaussian takes standard deviation rather than variance as its second argument.)

    W[0] = 0
    W[128] = np.random.gaussian(0, 128**0.5)

Then as described in the earlier post, we would generate W(64) by taking the average of W(0) and W(128) and adding a N(0, 32) sample. Note that the variance in the sample is 32, not 64. See the earlier post for an explanation.

    W[64] = 0.5*(W[0] + W[128]) + np.random.gaussian(0, 32**0.5)

Next we would fill in W(32) and W(96). To fill in W(32) we’d start by averaging W(0) and W(64), then add a N(0, 16) sample. And to fill in W(96) we’d average W(64) and W(128) and add a N(0, 16) sample.

    
    W[32] = 0.5*(W[0] + W[64] ) + np.random.gaussian(0, 16**0.5)
    W[64] = 0.5*(W[0] + W[128]) + np.random.gaussian(0, 16**0.5)

We would continue this process until we have samples at 0, 1, 2, …, 128. We won’t get the same Brownian path that we would have gotten if we’d sampled the paths in order—this is all random after all—but we’d get an instance of something created by the same stochastic process.

To state the fractal property explicitly, if W(t) is a Brownian path over [0, T], then

W(c² T)/c

is a Brownian path over [0, c² T]. Brownian motion looks the same at all scales.

***

The “less clear” part depends on who the reader is. The more concise version is more clear to someone accustomed to reading such code. But the loop more explicitly reflects the sequential nature of constructing the path.

Let W(t) be a standard Wiener process, a.k.a. a one-dimensional Brownian motion.

We can produce a discrete realization of W by first setting W(0) = 0. Then let W(1) be a sample from a N(0, 1) random variable. Then let W(2) be W(1) plus another N(0, 1) sample. At each integer n > 0, W(n) equals W(n-1) plus a N(0, 1) sample. Say we do this for n up to 1000.

Then for whatever reason we wish we’d produced a discretized realization with a sample at 800.5. We could start over, generating samples at steps of size 1/2. This time we’d add a N(0, 1/2) sample at each step [1].

We’d have to throw away a lot of work. How might we reuse the values we’ve already sampled?

(Update: Thanks to the comments I learned that interpolating Brownian motion is called a “Brownian bridge” in finance.)

The first thing you might think of might be to set W(800.5) to be the average of W(800) and W(801). That would be appropriate if we were sampling a deterministic function, but not when we’re sampling a stochastic process.

To decide the right thing to do we have to look back at the definition of a standard Wiener process. The three axioms are

1. W(0) = 0.
2. For s > t ≥ 0. W(s) – W(t) has distribution N(0, s – t).
3. For v ≥ u ≥ s ≥ t ≥ 0. W(s) – W(t) and W(v) – W(u) are independent.

Using linear interpolation violates (3) because if we set W(800.5) to be the average of W(800) and W(801), then W(801) – W(800.5) and W(800.5) – W(800) are not independent: they’re equal. That’s as dependent as it gets.

Another thing you might think of would be to add a N(0, 0.5) sample to W(800) because that’s exactly what you’d do if you had generated the Brownian path from scratch taking step sizes 0.5. And if you threw away the samples from 801 onward and started over again, that would work. But if you want to keep W(801) and its successors, you have to do something else.

Why wouldn’t this work? Because

W(801) – W(800.5)

would have too much variance, 1.5 when it should be 0.5.

The right thing to do is a sort of merger of the two ideas above: do the linear interpolation and add a N(0, 0.25) sample.

That is, we take

W(800.5) = 0.5 W(800) + 0.5 W(801) + N(0, 0.25).

Does this work? Let’s check that W(801) – W(800.5) has the right distribution.

W(801) – W(800.5) = 0.5W(801) – 0.5 W(800) – N(0, 0.25)

The distribution of

X = W(801) – W(800)

is N(0, 1), so 0.5 X has distribution N(0, 0.25). Thus

W(801) – W(800.5)

has the distribution of the difference of two N(0, 0.25) random variables, which is N(0, 0.5), as it should be.

A similar argument shows that W(800.5) – W(800) also has the right distribution.

Related posts

• Random walks and inheritance
• Random walk on quaternions
• Rapidly mixing random walks on graphs

[1] There’s always some ambiguity when you see a normal distribution with numeric parameters. Is the second parameter the standard deviation or the variance? If you see N(μ, σ) then by convention the second parameter σ is the standard deviation. If you see N(μ, σ²) then by convention the second parameter σ² is the variance. If you see N(0, 1/2) as above, is 1/2 a standard deviation or a variance? I intend it to be a variance.

This is an example of why some have called conventional statistical notation a “nomenclatural abomination.” You often can’t tell what a literal number means. You have to ask “If you were to write that number as a variable, what notation would you use?” In this case, you’d want to know whether 1/2 is a σ or a σ².

It’s well known that a binomial random variable can be approximated by a Poisson random variable, and under what circumstances the approximation is particularly good. See, for example, this post.

A binomial random variable is the sum of iid (independent, identically distributed) Bernoulli random variables. But what if the Bernoulli random variables don’t have the same distribution. That is, suppose you’re counting the number of heads seen in flipping n coins, where each coin has a potentially different probability of coming up heads. Will a Poisson approximation still work?

This post will cite three theorems on the error in approximating a sum of n independent Bernoulli random variables, each with a different probability of success p_i. I’ll state each theorem and very briefly discuss its advantages. The theorems can be found in [1].

Setup

For i = 1, 2, 3, …, n let X_i be Bernoulli random variables with

Prob(X_i = 1) = p_i

and let X with no subscript be their sum:

X = X₁ + X₂ + X₃ + … + X_n

We want to approximate the distribution of X with a Poisson distribution with parameter λ. We will measure the error in the Poisson approximation by the maximum difference between the mass density function for X and the mass density function for a Poisson(λ) random variable.

Sum of p‘s

We consider two ways to choose λ. The first is

λ = p₁ + p₂ + p₃ + … + p_n.

For this choice we have two different theorems that give upper bounds on the approximation error. One says that the error is bounded by the sum of the squares of the p‘s

p₁² + p₂² + p₃² + … + p_n²

and the other says it is bounded by 9 times the maximum of the p‘s

9 max(p₁, p₂, p₃, …,  p_n).

The sum of squares bound will be smaller when n is small and the maximum bound will be smaller when n is large.

Sum of transformed p‘s

The second way to choose λ is

λ = λ₁ + λ₂ + λ₃ + … + λ_n

where

λ_i = -log(1 – p_i).

In this case the bound on the error is one half the sum of the squared λ’s:

(λ₁² + λ₂² + λ₃² + … + λ_n²)/2.

When p_i is small, λ_i ≈ p_i. In this case the error bound for the transformed Poisson approximation will be about half that of the one above.

Related posts

• Normal approximation to binomial
• Camp-Paulson approximation to binomial
• Relative error in normal approximations

[1] R. J. Serfling. Some Elementary Results on Poisson Approximation in a Sequence of Bernoulli Trials. SIAM Review, Vol. 20, No. 3 (July, 1978), pp. 567-579.

It occurred to me recently that a problem I solved numerically years ago could be solved analytically, the problem of determining beta distribution parameters so that the distribution has a specified mean and variance. The calculation turns out to be fairly simple. Maybe someone has done it before.

Problem statement

The beta distribution has two positive parameters, a and b, and has probability density proportional to [1]

for x between 0 and 1.

The mean of a beta(a, b) distribution is

and the variance is

Given μ and σ² we want to solve for a and b. In order for the problem to be meaningful μ must be between 0 and 1, and σ² must be less than μ(1-μ). [2]

As we will see shortly, these two necessary conditions for a solution are also sufficient.

Visualization

Graphically, we want to find the intersection of a line of constant mean

with a line of constant variance.

Note that the scales in the two plots differ.

Solution

If we set

then b = ka and so we can eliminate b from the equation for variance to get

Now since a > 0, we can divide by a²

and from there solve for a:

and b = ka.

We require σ² to be less than μ(1-μ), or equivalently we require the ratio of μ(1-μ) to σ² to be greater than 1. It works out that the solution a is the product of the mean and the amount by which the ratio of μ(1-μ) to σ² exceeds 1.

Verification

Here is a little code to check for errors in the derivation above. It generates μ and σ² values at random, solves for a and b, then checks that the beta(a, b) distribution has the specified mean and variance.

    from scipy.stats import uniform, beta
    
    for _ in range(100):
        mu = uniform.rvs()
        sigma2 = uniform.rvs(0, mu*(1-mu))
    
        a = mu*(mu*(1-mu)/sigma2 - 1)
        b = a*(1-mu)/mu
    
        x = beta(a, b)
        assert( abs(x.mean() - mu) < 1e-10)
        assert( abs(x.var() - sigma2) < 1e-10)
    
    print("Done")

Related posts

• Determining a distribution from two quantiles
• Error in the normal approximation to a beta
• Diagram of probability distribution relationships

[1] It’s often easiest to think of probability densities ignoring proportionality constants. Densities integrate to 1, so the proportionality constants are determined by the rest of the expression for the density. In the case of the beta distribution, the proportionality constant works out to Γ(a + b) / Γ(a) Γ(b).

[2] The variance of a beta distribution factors into μ(1-μ)/(a + b + 1), so it is less than μ(1-μ).