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-μ).

Let Z be a standard normal random variable. Then there is a unique x such that

Pr(Z < x) = x.

That is, Φ has a unique fixed point where Φ is the CDF of a standard normal.

It’s easy to find the fixed point: start anywhere and iterate Φ.

Here’s a cobweb plot that shows how the iterates converge, starting with -2.

The black curve is a plot of Φ. The blue stairstep is a visualization of the iterates. The stairstep pattern comes from outputs turning into inputs. That is, vertical blue lines connect and input value x to its output value y, and the horizontal blue lines represent sliding a y value over to the dotted line y = x in order to turn it into the next x value.

y‘s become x‘s. The blue lines get short when we’re approaching the fixed point because now the outputs approximately equal the inputs.

Here’s a list of the first 10 iterates.

    0.02275
    0.50907
    0.69465
    0.75636
    0.77528
    0.78091
    0.78257
    0.78306
    0.78320
    0.78324

So it seems that after 10 iterations, we’ve converged to the fixed point in the first four decimal places.

A few days ago I wrote about the expected number of roots in a random polynomial where each coefficient is drawn from a standard normal, i.e. a Gaussian distribution with mean 0 and variance 1.

Another class of random polynomials, one that comes up in applications to physics, draws each coefficient from a different distribution. Specifically, the nth coefficient in an Nth degree polynomial is drawn from a normal distribution with mean zero and variance N choose n. As before, our reference is [1].

It turns out that for these random polynomials, the expected number of real zeros is √N. This is not an asymptotic approximation.

For an example, I created a random 44th degree polynomial.

    from scipy.stats import norm
    from scipy.special import binom

    N = 44
    coeff = [norm.rvs(scale=binom(N,n)**0.5) for n in range(N+1)]

I plotted the (tanh of) this polynomial as described here.

This polynomial clearly has 8 real roots. The expected number of roots for a 44th degree random polynomial of this type is √44 = 6.63.

[1] Alan Edelman and Eric Kostlan. How many zeros of a random polynomial are real? Bulletin of the AMS. Volume 32, Number 1, January 1995

A matrix is said to have a saddlepoint if an element is the smallest element in its row and the largest element in its column. For example, 0.546 is a saddlepoint in the matrix below because it is the smallest element in the third row and the largest element in the first column.

In [1], Goldman considers the probability of an m by n matrix having a saddlepoint if its elements are chosen independently at random from the same continuous distribution. He reports this probability as

but it may be easier to interpret as

When m and n are large, the binomial coefficient in the denominator is very large.

For square matrices, m and n are equal, and the probability above is 2n/(n + 1) divided by the nth Catalan number C_n.

A couple years ago I wrote about bounds on the middle binomial coefficient

by E. T. H. Wang. These bounds let us conclude that P(n, n), the probability of an n by n matrix having a saddlepoint, is bounded as follows.

Here’s a plot of P(n, n) and its bounds.

This shows that the probability goes to zero exponentially as n increases, being impossible to visually distinguish from 0 in the plot above when n = 8.

The upper bound on P(n, n) based on Wang’s inequalities is tighter than the corresponding lower bound. This is because Wang’s lower bound on the middle binomial coefficient is tighter (see plot here) and we’ve taken reciprocals.

Simulation

The theory above says that if we generate 4 by 3 matrices with independent normal random values, each matrix with have a saddle point with probability 1/5. The following code generates 10,000 random 4 by 3 matrices and counts what proportion have a saddle point.

    
    import numpy as np
    from scipy.stats import norm

    def has_saddlepoint(M):
        m, n = M.shape
        row_min_index = [np.argmin(M[i,:]) for i in range(m)]
        col_max_index = [np.argmax(M[:,j]) for j in range(n)]
        for i in range(m):
            if col_max_index[row_min_index[i]] == i:
                return True
        return False

    m, n = 4, 3
    N = 10000
    saddlepoints = 0

    for _ in range(N):
        M = norm.rvs(size=(m,n))
        if has_saddlepoint(M):
            saddlepoints += 1
    print(saddlepoints / N)

When I ran this, I got 0.1997, very close to the expected value.

More on random matrices

• Circular law for random matrices
• Eigenvalues of symmetric Gaussian matrices
• Spectra of random graphs

[1] A. J. Goldman. The Probability of a Saddlepoint. The American Mathematical Monthly, Vol. 64, No. 10, pp. 729-730

When the rule for stopping an experiment depends on the data in the experiment, the results could be biased if the stopping rule isn’t taken into account in the analysis [1].

For example, suppose Alice wants to convince Bob that π has a greater proportion of even digits than odd digits.

Alice: I’ll show you that π has more even digits than odd digits by looking at the first N digits. How big would you like N to be?

Bob: At least 1,000. Of course more data is always better.

Alice: Right. And how many more even than odd digits would you find convincing?

Bob: If there are at least 10 more evens than odds, I’ll believe you.

Alice: OK. If you look at the first 2589 digits, there are 13 more even digits than odd digits.

Now if Alice wanted to convince Bob that there are more odd digits, she could do that too. If you look at the first 2077 digits, 13 more are odd than even.

No matter what two numbers Bob gives, Alice can find a sample size that will give the result she wants. Here’s Alice’s Python code.

    from mpmath import mp
    import numpy as np

    N = 3000 
    mp.dps = N+2
    digits = str(mp.pi)[2:]

    parity = np.ones(N, dtype=int)
    for i in range(N):
        if digits[i] in ['1', '3', '5', '7', '9']:
            parity[i] = -1
    excess = parity.cumsum()
    print(excess[-1])
    print(np.where(excess == 13))
    print(np.where(excess == -13))

The number N is a guess at how far out she might have to look. If it doesn’t work, she increases it and runs the code again.

The array parity contains a 1 in positions where the digits of π (after the decimal point) are even and a -1 where they are odd. The cumulative sum shows how many more even than odd digits there have been up to a given point, a negative number meaning there have been more odd digits.

Alice thought that stopping when there are exactly 10 more of the parity she wants would look suspicious, so she looked for places where the difference was 13.

Here are the results:

    [ 126,  128,  134,  …,  536, 2588, … 2726]
    [ 772,  778,  780,  …,  886, 2076, … 2994]

There’s one minor gotcha. The array excess is indexed from zero, so Alice reports 2589 rather than 2588 because the 2589th digit has index 2588.

Bob’s mistake was that he specified a minimum sample size. By saying “at least 1,000” he gave Alice the freedom to pick the sample size to get the result she wanted. If he specified an exact sample size, there probably would be either more even digits or more odd digits, but there couldn’t be both. And if he were more sophisticated enough, he could pick an excess value that would be unlikely given that sample size.

Related posts

• Stopping trials of ineffective drugs sooner
• Finding coffee in pi
• Balancing profit and learning in A/B testing

[1] This does not contradict the likelihood principle; it says that informative stopping rules should be incorporated into the likelihood function.

The Fibonacci numbers are defined by F₁ = F₂ = 1, and for n > 2,

F_n = F_n-1 + F_n-2.

A random Fibonacci sequence f is defined similarly, except the addition above is replaced with a subtraction with probability 1/2. That is, f₁ = f₂ = 1, and for n > 2,

f_n = f_n-1 + b f_n-2

where b is +1 or -1, each with equal probability.

Here’s a graph a three random Fibonacci sequences.

Here’s the Python code that was used to produce the sequences above.

    import numpy as np

    def rand_fib(length):
        f = np.ones(length)
        for i in range(2, length):
            b = np.random.choice((-1,1))
            f[i] = f[i-1] + b*f[i-2]
        return f

It’s easy to see that the nth random Fibonacci number can be as large as the nth ordinary Fibonacci number if all the signs happen to be positive. But the numbers are typically much smaller.

The nth (ordinary) Fibonacci number asymptotically approaches φⁿ is the golden ratio, φ = (1 + √5)/2 = 1.618…

Another way to say this is that

The nth random Fibonacci number does not have an asymptotic value—it wanders randomly between positive and negative values—but with probability 1, the absolute values satisfy

This was proved in 1960 [1].

Here’s a little Python code to show that we get results consistent with this result using simulation.

    N = 500
    x = [abs(rand_fib(N)[-1])**(1/N) for _ in range(10)]
    print(f"{np.mean(x)} ± {np.std(x)}")

This produced

1.1323 ± 0.0192

which includes the theoretical value 1.1320.

Update: The next post looks at whether every integer appear in a random Fibonacci sequence. Empirical evidence suggests the answer is yes.

Related posts

• Fibonacci meets Pythagoras
• Fibonacci meets Chebyshev
• Power method and Fibonacci

[1] Furstenberg and Kesten. Products of random matrices. Ann. Math. Stat. 31, 457-469.