Uniform distribution of powers mod 1

A few days ago I wrote about how powers of the golden ratio are nearly integers but powers of π are not. This post is similar but takes a little different perspective. Instead of looking at how close powers are to the nearest integers, we’ll look at how close they are to their floor, the largest integer below. Put another way, we’ll throw away the integer parts and look at the decimal parts.

First a theorem:

For almost all x > 1, the sequence (xn) for n = 1, 2, 3, … is u.d. mod 1. [1]

Here “almost all” is a technical term meaning that the set of x‘s for which the statement above does not hold has Lebesgue measure zero. The abbreviation “u.d.” stands for “uniformly distributed.” A sequence uniformly distributed mod 1 if the fractional parts of the sequence are distributed like uniform random variables.

Even though the statement holds for almost all x, it’s hard to prove for particular values of x. And it’s easy to find particular values of x for which the theorem does not hold.

From [1]:

… it is interesting to note that one does not know whether sequences such as (en), (πn), or even ((3/2)n) are u.d. mod 1 or not.

Obviously powers of integers are not u.d. mod 1 because their fractional parts are all 0. And we’ve shown before that powers of the golden ratio and powers of the plastic constant are near integers, i.e. their fractional parts cluster near 0 and 1.

The curious part about the quote above is that it’s not clear whether powers of 3/2 are uniformly distributed mod 1. I wouldn’t expect powers of any rational number to be u.d. mod 1. Either my intuition was wrong, or it’s right but hasn’t been proved, at least not when [1] was written.

The next post will look at powers of 3/2 mod 1 and whether they appear to be uniformly distributed.

* * *

[1] Kuipers and Niederreiter, Uniform Distribution of Sequences

Plastic powers

Last week I wrote a blog post showing that powers of the golden ratio are nearly integers. Specifically, the distance from φn to the nearest integer decreases exponentially as n increases. Several people pointed out that the golden constant is a Pisot number, the general class of numbers whose powers are exponentially close to integers.

The so-called plastic constant P is another Pisot number, in fact the smallest Pisot number. P is the real root of x3x – 1 = 0.

P = \frac{ (9 - \sqrt{69})^{1/3} + (9 + \sqrt{69})^{1/3} }{ 2^{1/3} \,\,\, 3^{2/3} } = 1.3247\ldots

Because P is a Pisot number, we know that its powers will be close to integers, just like powers of the golden ratio, but the way they approach integers is more interesting. The convergence is slower and less regular.

We will the first few powers of P, first looking at the distance to the nearest integer on a linear scale, then looking at the absolute value of the distance on a logarithmic scale.

distance from powers of plastic constant to nearest integer

distance from powers of plastic constant to nearest integer, log scale

As a reminder, here’s what the corresponding plots looked like for the golden ratio.

distance from powers of golden ratio to nearest integer

distance from powers of golden ratio to nearest integer, log scale

Golden powers are nearly integers

Nautilus, golden ratio

This morning I was reading Terry Tao’s overview of the work of Yves Meyer and ran across this line:

The powers φ, φ2, φ3, … of the golden ratio lie unexpectedly close to integers: for instance, φ11 = 199.005… is unusually close to 199.

I’d never heard that before, so I wrote a little code to see just how close golden powers are to integers.

Here’s a plot of the difference between φn and the nearest integer:

distance from powers of golden ratio to nearest integer

(Note that if you want to try this yourself, you need extended precision. Otherwise you’ll get strange numerical artifacts once φn is too large to represent exactly.)

By contrast, if we make the analogous plot replacing φ with π we see that the distance to the nearest integer looks like a uniform random variable:

distance from powers of pi to nearest integer

The distance from powers of φ to the nearest integer decreases so fast that cannot see it in the graph for moderate sized n, which suggests we plot the difference on the log scale. (In fact we plot the log of the absolute value of the difference since the difference could be negative and the log undefined.) Here’s what we get:

absolute distance from powers of golden ratio to nearest integer on log scale

After an initial rise, the curve is apparently a straight line on a log scale, i.e. the absolute distance to the nearest integer decreases almost exactly exponentially.

Related posts:

Inverse Fibonacci numbers

As with the previous post, this post is a spinoff of a blog post by Brian Hayes. He considers the problem of determining whether a number n is a Fibonacci number and links to a paper by Gessel that gives a very simple solution: A positive integer n is a Fibonacci number if and only if either 5n2 – 4 or 5n2 + 4 is a perfect square.

If we know n is a Fibonacci number, how can we tell which one it is? That is, if n = Fm, how can we find m?

For large m, Fm is approximately φm / √ 5 and the error decreases exponentially with m. By taking logs, we can solve for m and round the result to the nearest integer.

We can illustrate this with SymPy. First, let’s get a Fibonacci number.

      >>> from sympy import *
      >>> F = fibonacci(100)
      >>> F
      354224848179261915075

Now let’s forget that we know F is a Fibonacci number and test whether it is one.

      >>> sqrt(5*F**2 - 4)
      sqrt(627376215338105766356982006981782561278121)

Apparently 5F2 – 4 is not a perfect square. Now let’s try 5F2 + 4.

      >>> sqrt(5*F**2 + 4)
      792070839848372253127

Bingo! Now that we know it’s a Fibonacci number, which one is it?

      >>> N((0.5*log(5) + log(F))/log(GoldenRatio), 10)
      100.0000000

So F must be the 100th Fibonacci number.

It looks like our approximation gave an exact result, but if we ask for more precision we see that it did not.

      >>> N((0.5*log(5) + log(F))/log(GoldenRatio), 50)
      99.999999999999999999999999999999999999999996687654

Related posts:

Infinite primes via Fibonacci numbers

A well-known result about Fibonacci numbers says

gcd(Fm, Fn) = Fgcd(m, n)

In words, the greatest common divisor of the mth and nth Fibonacci numbers is the gth Fibonacci number where g is the greatest common divisor of m and n. You can find a proof here.

M. Wunderlich used this identity to create a short, clever proof that there are infinitely many primes.

Suppose on the contrary that there are only k prime numbers, and consider the set of Fibonacci numbers with prime indices: Fp1, Fp2, … Fpk. The Fibonacci theorem above tells us that these Fibonacci numbers are pairwise relatively prime: each pair has greatest common divisor F1 = 1.

Each of the Fibonacci numbers in our list must have only one prime factor. Why? Because we have assumed there are only k prime numbers, and no two of our Fibonacci numbers share a prime factor. But F19 = 4181 = 37*113. We’ve reached a contradiction, and so there must be infinitely many primes.

Source: M. Wunderlich, Another proof of the infinite primes theorem. American Mathematical Monthly, Vol. 72, No. 3 (March 1965), p. 305.

Here are a couple more unusual proofs that there are infinitely many primes. The first uses a product of sines. The second is from Paul Erdős. It also gives a lower bound on π(N), the number of primes less than N.

Three proofs that 2017 is prime

Aaron Meurer asked on Twitter whether there’s a proof that 2017 is prime that would fit into 140 characters.

My first reply was this:

sqrt(2017) < 45.
2017 not divisible by 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, or 43.
Ergo prime.

I’m not sure that’s what he had in mind. There’s some implied calculation (which I didn’t actually do), so it’s kinda cheating. It would be interesting if there were something special about 2017 that would allow a more transparent proof.

(Implicit in the proof above is the theorem that if a number has a prime factor, it has a prime factor less than it’s square root. If you multiply together numbers bigger than the square root of n, you get a product bigger than n.)

Then for fun I gave two more proofs that are more sophisticated but would require far too much work to actually carry out by hand.

The first uses Fermat’s little theorem:

For 0 < a < 2017, a2016 – 1 is divisible by 2017.
2017 is not one of the three Carmichael numbers < 2465.
Ergo prime.

Fermat’s little theorem says that if p is prime, then for any 0 < ap, ap – 1 – 1 is divisible by p. This is actually an efficient way to prove that a number is not prime. Find a number a such that the result doesn’t hold, and you’ve proved that p isn’t prime. For small numbers, the easiest way to show a number is not prime is to show its factors. But for very large numbers, such as those used in cryptography, it’s efficient to have a way to prove that a number has factors without having to actually produce those factors.

Unfortunately, Fermat’s little theorem gives a necessary condition for a number to be prime, but not a sufficient condition. It can appear to be prime for every witness (the bases a are called witnesses) and still not be a prime. The Carmichael numbers pass the Fermat primailty test without being prime. The first four are 561, 1105, 1729, and 2465.

For more on using Fermat’s little theorem to test for primality, see Probability that a number is prime.

The final proof was this:

2016! + 1 is divisible by 2017, and so by Wilson’s theorem 2017 is prime.

Unlike Fermat’s little theorem, Wilson’s theorem gives necessary and sufficient conditions for a number to be prime. A number n is prime if and only if (n-1)! + 1 is divisible by n. In theory you could use Wilson’s theorem to test whether a number is prime, but this would be horrendously inefficient. 2016! has 5,789 digits. (You can find out how many digits n! has without computing it using a trick described here.)

Despite its inefficiency, you can actually use Wilson’s theorem and SymPy to prove that 2017 is prime.

      >>> from sympy import factorial
      >>> x = factorial(2016) + 1
      >>> x % 2017
      0

 

A short, unusual proof that there are infinitely many primes

Sam Northshield [1] came up with the following clever proof that there are infinitely many primes.

Suppose there are only finitely many primes and let P be their product. Then

0 < \prod_p \sin\left( \frac{\pi}{p} \right) = \prod_p \sin\left(\frac{\pi(1+2P)}{p} \right) = 0

The original publication gives the calculation above with no explanation. Here’s a little commentary to explain the calculation.

Since prime numbers are greater than 1, sin(π/p) is positive for every prime. And a finite product of positive terms is positive. (An infinite product of positive terms could converge to zero.)

Since p is a factor of P, the arguments of sine in the second product differ from those in the first product by an integer multiple of 2π, so the corresponding terms in the two products are the same.

There must be some p that divides 1 + 2P, and that value of p contributes the sine of an integer multiple of π to the product, i.e. a zero. Since one of the terms in the product is zero, the product is zero. And since zero is not greater than zero, we have a contradiction.

* * *

[1] A One-Line Proof of the Infinitude of Primes, The American Mathematical Monthly, Vol. 122, No. 5 (May 2015), p. 466

Computing discrete logarithms with baby-step giant-step algorithm

At first “discrete logarithm” sounds like a contradiction in terms. Logarithms aren’t discrete, not as we usually think of them. But you can define and compute logarithms in modular arithmetic.

What is a logarithm? It’s the solution to an exponential equation. For example, the logarithm base 10 of 2 is the solution to the equation 10x = 2, and so x =0.30103. Similarly, you could look for the logarithm base 10 of 2 modulo 19. This is an integer value of x such that 10x = 2 mod 19, and you can show that 17 is a solution.

If we work modulo an integer n, the discrete logarithm base a of a number y is an integer x such that ax = y. For some values of a there may not be a solution, but there will be a solution if a is a generator of the integers mod n.

Brute force the simplest algorithm for finding discrete logarithms: try x = 1, 2, 3, …, n until you find a value of x satisfying ax = y. The problem with brute force is that the expected run time is on the order of n, and n is often very large in application.

Discrete logarithms are expensive to compute, but we can do better than brute force. (Cryptographic applications take advantage of the fact that discrete logarithms are hard to compute.) There’s a simple algorithm by Daniel Shanks, known as the baby-step giant-step algorithm, that reduces the run time from order n to order roughly √n. (Actually O(√n log n) for reasons we’ll see soon.)

Let s be the ceiling of the square root of n. The idea is to search for x by taking giant steps forward (multiples of s) and baby (single) steps backward.

Let G be the set of pairs (ags mod n, gs) where g ranges from 1 to s. These are the giant steps.

Let B be the set of pairs (yab mod n, b) where b ranges from 0 to s-1. These are the baby steps.

Now look for a pair in G and a pair in B with the matching first components, i.e. yab = ags. Then x = gsb is the required solution.

Constructing the sets G and requires O(s) operations, and matching the two sets takes O(s log s) operations.

Here’s an example, going back to the problem above of finding the discrete log base 10 of 2 mod 19, using a little Python code to help with the calculations.

The square root of 19 is between 4 and 5 so s = 5.

     >>> [(2*10**r % 19, r) for r in range(5)]
     [(2, 0), (1, 1), (10, 2), (5, 3), (12, 4)]
     >>> [(10**(4*t) % 19, 4*t) for t in range(1,6)]
     [(6, 4), (17, 8), (7, 12), (4, 16), (5, 20)]

The number 5 appears as the first element of a pair in both B and G. We have (5, 3) in B and (5, 20) in G so x = 20 – 3 = 17.

Related: Applied number theory

Periods of fractions

Suppose you have a fraction a/b where 0 < ab, and a and b are relatively prime integers. The decimal expansion of a/b either terminates or it has an initial non-repeating part followed by a repeating part.

How long is the non-repeating part? How long is the period of the repeating part?

The answer depends on the prime factorization of the denominator b. If b has the form

b = 2α 5β

then the decimal expansion has r digits where r = max(α, β).

Otherwise b has the factorization

b = 2α 5β d

where d is some integer greater than 1 and relatively prime to 10. Then the decimal expansion of a/b has r non-repeating digits and a repeating part with period s where r is as above, and s is the smallest positive integer such that

d | 10s– 1,

i.e. the smallest s such that 10s – 1 is divisible by d. How do we know there exists any such integer s? This isn’t obvious.

Fermat’s little theorem tells us that

d | 10d – 1 – 1

and so we could take s = d – 1, though this may not be the smallest such s. So not only does s exist, we know that it is at most d – 1. This means that the period of the repeating part of a/b is no more than d – 1 where d is what’s left after factoring out as many 2’s and 5’s from b as possible.

By the way, this means that you can take any integer d, not divisible by 2 or 5, and find some integer k such that dk consists only of 9’s.

Related post: Cyclic fractions

Münchausen numbers

Baron Münchausen

Baron Münchausen

The number 3435 has the following curious property:

3435 = 33 + 44 + 33 + 55.

It is called a Münchausen number, an allusion to fictional Baron Münchausen. When each digit is raised to its own power and summed, you get the original number back. The only other Münchausen number is 1.

At least in base 10. You could look at Münchausen numbers in other bases. If you write out a number n in base b, raise each of its “digits” to its own power, take the sum, and get n back, you have a Münchausen number in base b. For example 28 is a Münchausen number in base 9 because

28ten = 31nine = 33 + 11

Daan van Berkel proved that there are only finitely many Münchausen in any given base. In fact, he shows that a Münchausen number in base b cannot be greater than 2bb, and so you could do a brute-force search to find all the Münchausen numbers in any base.

The upper bound 2bb grows very quickly with b and so brute force becomes impractical for large b. If you wanted to find all the hexadecimal Münchausen numbers you’d have to search 2*1616 = 36,893,488,147,419,103,232 numbers. How could you do this more efficiently?

Prime factors, phone numbers, and the normal distribution

Telephone numbers typically have around three distinct prime factors.

The length of a phone number varies by country, but US a phone number is a 10 digit number, and 10-digit numbers are often divisible by three different prime numbers, give or take a couple. Assuming phone numbers are scattered among possible 10-digit numbers in a way that doesn’t bias their number of prime factors, these numbers will often have between 1 and 5 prime factors. If a country has 9- or 11-digit phone numbers, the result is essentially the same.

Let ω(n) be the number of distinct prime factors of n. Then the Erdős–Kac theorem says roughly that ω(n) is distributed like a normal random variable with mean and variance log log n. More precisely, fix two numbers a and b.  For a given value of x, count the proportion of positive integers less than x where

(ω(n) – log log n) / sqrt( log log n)

is between a and b. Then in the limit as x goes to infinity, this proportion approaches the probability that a standard normal random variable is between a and b.

So by that heuristic, the number of distinct prime factors of a 10-digit number is approximately normally distributed with mean and variance log log 10^11 = 3.232, and such a distribution is between 1 and 5 around 73% of the time.

My business phone number, for example, is 8324228646. Obviously this is divisible by 2. In fact it equals 2 × 32 × 462457147, and so it has exactly three distinct prime factors: 2, 3, and 462457147.

Here’s how you could play with this using Python.

    from sympy.ntheory import factorint

    def omega(n):
        return len(factorint(n))

I looked in SymPy and didn’t see an implementation of ω(n) directly, but it does have a function factorint that returns the prime factors of a number, along with their multiplicities, in a dictionary. So ω(n) is just the size of that dictionary.

I took the first 20 phone numbers in my contacts and ran them through omega and got results consistent with what you’d expect from the theory above. One was prime, and none had more than five factors.

Bar chart of umber of prime factors in a sample of phone numbers with heights [1, 3, 5, 8, 3]

Benford’s law, chi-square, and factorials

A while back I wrote about how the leading digits of factorials follow Benford’s law. That is, if you look at just the first digit of a sequence of factorials, they are not evenly distributed. Instead, 1’s are most popular, then 2’s, etc. Specifically, the proportion of factorials starting with n is roughly log10(1 + 1/n).

Someone has proved that the limiting distribution of leading digits of factorials exactly satisfies Benford’s law. But if we didn’t know this, we might use a chi-square statistic to measure how well the empirical results match expectations. As I argued in the previous post, statistical tests don’t apply here, but they can be useful anyway in giving us a way to measure the size of the deviations from theory.

Benford’s law makes a better illustration of the chi-square test than the example of prime remainders because the bins are unevenly sized, which they’re allowed to be in general. In the prime remainder post, they were all the same size.

The original post on leading digits of factorial explains why we compute the leading digits the way we do. Only one detail has changed: the original post used Python 2 and this one uses Python 3. Integer division was the default in Python 2, but now in Python 3 we have to use // to explicitly ask for integer division, floating point division being the new default.

Here’s a plot of the distribution of the leading digits for the first 500 factorials.

And here’s code to compute the chi-square statistic:

    from math import factorial, log10

    def leading_digit_int(n):
        while n > 9:
            n = n//10
        return n

    def chisq_stat(O, E):
        return sum( [(o - e)**2/e for (o, e) in zip(O, E)] )

    # Waste the 0th slot to make the code simpler.
    digits = [0]*10

    N = 500
    for i in range(N):
        digits[ leading_digit_int( factorial(i) ) ] += 1

    expected = [ N*log10(1 + 1/n) for n in range(1, 10) ]

    print( chisq_stat(digits[1:], expected) )

This gives a chi-square statistic of 7.693, very near the mean value of 8 for a chi-square distribution with eight degrees of freedom. (There are eight degrees of freedom, not nine, because if we know how many factorials start with the digits 1 through 8, we know how many start with 9.)

So the chi-square statistic suggests that the deviation from Benford’s law is just what we’d expect from random data following Benford’s law. And as we said before, this suggestion turns out to be correct.

Related posts:

Hypothesis testing and number theory

This post uses a hypothesis test for proportions to look at a couple conjectures in number theory. It is similar to my earlier post on the chi-square test and prime remainders. You could read this as a post on statistics or a post on number theory, depending on which you’re less familiar with.

Using statistical tests on number theory problems is kind of odd. There’s nothing random going on, so in that since the whole enterprise is unjustified. Nevertheless, statistical tests can be suggestive. They certainly don’t prove theorems, but they can give reason to suspect a theorem is true or false. In that sense, applying statistical tests to number theory isn’t all that different from applying them to more traditional settings.

First we’ll look at the remainders of primes modulo 4. Except for 2, all primes are odd, and so they either have remainder 1 or 3 when divided by 4. Brian Hayes wrote recently that Chebyshev noticed in the 1850’s that there seems to be more primes with remainder 3. Is the imbalance larger than one would expect to see from fair coin tosses?

Here’s some Python code to find the proportion of the first million primes (after 2) that have remainder 3 when divided by 4.

    from sympy import prime
    from scipy import sqrt

    n = 1000000
    rem3 = 0
    for i in range (2, n+2):
        if prime(i) % 4 == 3:
            rem3 += 1
    p_hat = rem3/n

This shows that of the first million odd primes, 500,202 are congruent to 3 mod 4. Would it be unusual for a coin to come up heads this many times in a million flips? To find out we’d compute the z-statistic:

z = \frac{\hat{p} - p}{\sqrt{pq/n}}

Here p is the proportion under the null hypothesis, q = 1 – p, and n is the sample size. In our case, the null hypothesis is p = 0.5 and n = 1,000,000. [1]

The code

    p = 0.5
    q = 1 - p
    z = (p_hat - p)/sqrt(p*q/n)

shows that z = 0.404, hardly a suspiciously large value. If we were looking at random values we’d see a z-score this large or larger 34% of the time. So in this case doesn’t suggest much.

Click to learn more about Bayesian statistics consulting

 
[1] The derivation of the z statistic is fairly quick. If the proportion of successes is p, then the number of successes out of n trials is binomial(np). For large n, this is has approximately the same distribution as a normal distribution with the same mean and variance, mean np and variance npq. The proportion of successes then has approximately mean p and standard deviation √(pq/n). Subtracting the mean and dividing by the standard deviation normalizes the distribution to have mean 0 and variance 1. So under the null hypothesis, the z statistic has a standard normal distribution.

Prime remainders too evenly distributed

First Brian Hayes wrote an excellent post about the remainders when primes are divided by other primes. Then I wrote a follow-on just focusing on the first part of his post. He mostly looked at pairs of primes, but I wanted to look in more detail at the first part of his post, simulating dice rolls by keeping the remainder when consecutive primes are divided by a fixed prime. For example, using a sequence of primes larger than 7 and taking their remainder by 7 to create values from 1 to 6.

The results are evenly distributed with some variation, just like dice rolls. In fact, given these results and results from a set of real dice rolls, most people would probably think the former are real because they’re more evenly distributed. A chi-squared goodness of fit test shows that the results are too evenly distributed compared to real dice rolls.

At the end of my previous post, I very briefly discuss what happens when you look a “dice” with more than six sides. Here I’ll go into a little more detail and look at a large number of examples.

In short, you either get a suspiciously good fit or a terrible fit. If you look at the remainder when dividing primes by m, you get values between 1 and m-1. You can’t get a remainder of 0 because primes aren’t divisible by m (or anything else!). If m itself is prime, then you get all the numbers between 1 and m-1, and as we’ll show below you get them in very even proportions. But if m isn’t prime, there are some remainders you can’t get.

The sequence of remainders looks random in the sense of being unpredictable. (Of course it is predictable by the algorithm that generates them, but it’s not predictable in the sense that you could look at the sequence out of context and guess what’s coming next.) The sequence is biased, and that’s the big news. Pairs of consecutive primes have correlated remainders. But I’m just interested in showing a different departure from a uniform distribution, namely that the results are too evenly distributed compared to random sequences.

The table below gives the chi-square statistic and p-value for each of several primes. For each prime p, we take remainders mod p of the next million primes after p and compute the chi-square goodness of fit statistic with p-2 degrees of freedom. (Why p-2? There are p-1 different remainders, and the chi-square test for k possibilities has k-1 degrees of freedom.)

The p-value column gives the probability of seeing at fit this good or better from uniform random data. (The p in p-value is unrelated to our use of p to denote a prime. It’s an unfortunate convention of statistics that everything is denoted p.) After the first few primes, the p-values are extremely small, indicating that such an even distribution of values would be astonishing from random data.

|-------+------------+------------|
| Prime | Chi-square | p-value    |
|-------+------------+------------|
|     3 |     0.0585 |   2.88e-02 |
|     5 |     0.0660 |   5.32e-04 |
|     7 |     0.0186 |   1.32e-07 |
|    11 |     0.2468 |   2.15e-07 |
|    13 |     0.3934 |   6.79e-08 |
|    17 |     0.5633 |   7.64e-10 |
|    19 |     1.3127 |   3.45e-08 |
|    23 |     1.1351 |   2.93e-11 |
|    29 |     1.9740 |   3.80e-12 |
|    31 |     2.0052 |   3.11e-13 |
|    37 |     2.5586 |   3.92e-15 |
|    41 |     3.1821 |   9.78e-16 |
|    43 |     4.4765 |   5.17e-14 |
|    47 |     3.7142 |   9.97e-18 |
|    53 |     3.7043 |   3.80e-21 |
|    59 |     7.0134 |   2.43e-17 |
|    61 |     5.1461 |   6.45e-22 |
|    67 |     7.1037 |   5.38e-21 |
|    71 |     7.6626 |   6.13e-22 |
|    73 |     7.5545 |   4.11e-23 |
|    79 |     8.0275 |   3.40e-25 |
|    83 |    12.1233 |   9.92e-21 |
|    89 |    11.4111 |   2.71e-24 |
|    97 |    12.4057 |   2.06e-26 |
|   101 |    11.8201 |   3.82e-29 |
|   103 |    14.4733 |   3.69e-26 |
|   107 |    13.8520 |   9.24e-29 |
|   109 |    16.7674 |   8.56e-26 |
|   113 |    15.0897 |   1.20e-29 |
|   127 |    16.4376 |   6.69e-34 |
|   131 |    19.2023 |   6.80e-32 |
|   137 |    19.1728 |   1.81e-34 |
|   139 |    22.2992 |   1.82e-31 |
|   149 |    22.8107 |   6.67e-35 |
|   151 |    22.8993 |   1.29e-35 |
|   157 |    30.1726 |   2.60e-30 |
|   163 |    26.5702 |   3.43e-36 |
|   167 |    28.9628 |   3.49e-35 |
|   173 |    31.5647 |   7.78e-35 |
|   179 |    33.3494 |   2.46e-35 |
|   181 |    36.3610 |   2.47e-33 |
|   191 |    29.1131 |   1.68e-44 |
|   193 |    29.9492 |   2.55e-44 |
|   197 |    34.2279 |   3.49e-41 |
|   199 |    36.7055 |   1.79e-39 |
|   211 |    41.0392 |   8.42e-40 |
|   223 |    39.6699 |   1.73e-45 |
|   227 |    42.3420 |   2.26e-44 |
|   229 |    37.1896 |   2.02e-50 |
|   233 |    45.0111 |   4.50e-44 |
|   239 |    43.8145 |   2.27e-47 |
|   241 |    51.3011 |   1.69e-41 |
|   251 |    47.8670 |   6.28e-48 |
|   257 |    44.4022 |   1.54e-53 |
|   263 |    51.5905 |   7.50e-49 |
|   269 |    59.8398 |   3.92e-44 |
|   271 |    59.6326 |   6.02e-45 |
|   277 |    52.2383 |   2.80e-53 |
|   281 |    52.4748 |   1.63e-54 |
|   283 |    64.4001 |   2.86e-45 |
|   293 |    59.7095 |   2.59e-52 |
|   307 |    65.2644 |   1.64e-52 |
|   311 |    63.1488 |   1.26e-55 |
|   313 |    68.6085 |   7.07e-52 |
|   317 |    63.4099 |   1.72e-57 |
|   331 |    66.3142 |   7.20e-60 |
|   337 |    70.2918 |   1.38e-58 |
|   347 |    71.3334 |   3.83e-61 |
|   349 |    75.8101 |   3.38e-58 |
|   353 |    74.7747 |   2.33e-60 |
|   359 |    80.8957 |   1.35e-57 |
|   367 |    88.7827 |   1.63e-54 |
|   373 |    92.5027 |   7.32e-54 |
|   379 |    86.4056 |   5.67e-60 |
|   383 |    74.2349 |   3.13e-71 |
|   389 |   101.7328 |   9.20e-53 |
|   397 |    86.9403 |   1.96e-65 |
|   401 |    90.3736 |   3.90e-64 |
|   409 |    92.3426 |   2.93e-65 |
|   419 |    95.9756 |   8.42e-66 |
|   421 |    91.1197 |   3.95e-70 |
|   431 |   100.3389 |   1.79e-66 |
|   433 |    95.7909 |   1.77e-70 |
|   439 |    96.2274 |   4.09e-72 |
|   443 |   103.6848 |   6.96e-68 |
|   449 |   105.2126 |   1.07e-68 |
|   457 |   111.9310 |   1.49e-66 |
|   461 |   106.1544 |   7.96e-72 |
|   463 |   116.3193 |   1.74e-65 |
|   467 |   116.2824 |   1.02e-66 |
|   479 |   104.2246 |   3.92e-79 |
|   487 |   116.4034 |   9.12e-73 |
|   491 |   127.2121 |   6.69e-67 |
|   499 |   130.9234 |   5.90e-67 |
|   503 |   118.4955 |   2.60e-76 |
|   509 |   130.9212 |   6.91e-70 |
|   521 |   118.6699 |   6.61e-82 |
|   523 |   135.4400 |   3.43e-71 |
|   541 |   135.9210 |   3.13e-76 |
|   547 |   120.0327 |   2.41e-89 |
|-------+------------+------------|

Chi-square goodness of fit test example with primes

chi squared

Yesterday Brian Hayes wrote a post about the distribution of primes. He showed how you could take the remainder when primes are divided by 7 and produce something that looks like rolls of six-sided dice. Here we apply the chi-square goodness of fit test to show that the rolls are too evenly distributed to mimic randomness. This post does not assume you’ve seen the chi-square test before, so it serves as an introduction to this goodness of fit test.

In Brian Hayes’ post, he looks at the remainder when consecutive primes are divided by 7, starting with 11. Why 11? Because it’s the smallest prime bigger than 7. Since no prime is divisible by any other prime, all the primes after 7 will have a remainder of between 1 and 6 inclusive when divided by 7. So the results are analogous to rolling six-sided dice.

The following Python code looks at prime remainders and (pseudo)random rolls of dice and computes the chi-square statistic for both.

First, we import some functions we’ll need.

    from sympy import prime
    from random import random
    from math import ceil

The function prime takes an argument n and returns the nth prime. The function random produces a pseudorandom number between 0 and 1. The ceiling function ceil rounds its argument up to an integer. We’ll use it to convert the output of random into dice rolls.

In this example we’ll use six-sided dice, but you could change num_sides to simulate other kinds of dice. With six-sided dice, we divide by 7, and we start our primes with the fifth prime, 11.

    num_sides = 6
    modulus = num_sides + 1

    # Find the index of the smallest prime bigger than num_sides
    index = 1
    while prime(index) <= modulus:
        index += 1

We’re going to take a million samples and count how many times we see 1, 2, …, 6. We’ll keep track of our results in an array of length 7, wasting a little bit of space since the 0th slot will always be 0. (Because the remainder when dividing a prime by a smaller number is always positive.)

    # Number of samples
    N = 1000000
    
    observed_primes = [0]*modulus
    observed_random = [0]*modulus

Next we “roll” our dice two ways, using prime remainders and using a pseudorandom number generator.

    for i in range(index, N+index):
        m = prime(i) % modulus
        observed_primes[m] += 1
        m = int(ceil(random()*num_sides))
        observed_random[m] += 1

The chi-square goodness of fit test depends on the observed number of events in each cell and the expected number. We expect 1/6th of the rolls to land in cell 1, 2, …, 6 for both the primes and the random numbers. But in a general application of the chi-square test, you could have a different expected number of observations in each cell.

    expected = [N/num_sides for i in range(1, modulus)]

The chi-square test statistic sums (O – E)2/E over all cells, where O stands for “observed” and E stands for “expected.”

    def chisq_stat(O, E):
        return sum( [(o - e)**2/e for (o, e) in zip(O, E)] )

Finally, we compute the chi-square statistic for both methods.

    ch = chisq_stat(observed_primes[1:], expected[1:])
    print(ch)

    ch = chisq_stat(observed_random[1:], expected[1:])
    print(ch)

Note that we chop off the first element of the observed and expected lists to get rid of the 0th element that we didn’t use.

When I ran this I got 0.01865 for the prime method and 5.0243 for the random method. Your results for the prime method should be the same, though you might have a different result for the random method.

Now, how do we interpret these results? Since we have six possible outcomes, our test statistics has a chi-square distribution with five degrees of freedom. It’s one less than the number of possibilities because the total counts have to sum to N; if you know how many times 1, 2, 3, 4, and 5 came up, you can calculate how many times 6 came up.

A chi-square distribution with ν degrees of freedom has expected value ν. In our case, we expect a value around 5, and so the chi-square value of 5.0243 is unremarkable. But the value of 0.01864 is remarkably small. A large chi-square statistics would indicate a poor fit, the observed numbers being suspiciously far from their expected values. But a small chi-square value suggests the fit is suspiciously good, closer to the expected values than we’d expect of a random process.

We can be precise about how common or unusual a chi-square statistic is by computing the probability that a sample from the chi square distribution would be larger or smaller. The cdf gives the probability of seeing a value this small or smaller, i.e. a fit this good or better. The sf gives the probability of seeing a value this larger or larger, i.e. a fit this bad or worse. (The scipy library uses sf for “survival function,” another name for the ccdf, complementary cumulative distribution function).

    from scipy.stats import chi2
    print(chi2.cdf(ch, num_sides-1), chi2.sf(ch, num_sides-1))

This says that for the random rolls, there’s about a 41% chance of seeing a better fit and a 59% chance of seeing a worse fit. Unremarkable.

But it says there’s only a 2.5 in a million chance of seeing a better fit than we get with prime numbers. The fit is suspiciously good. In a sense this is not surprising: prime numbers are not random! And yet in another sense it is surprising since there’s a heuristic that says primes act like random numbers unless there’s a good reason why in some context they don’t. This departure from randomness is the subject of research published just this year.

If you look at dice with 4 or 12 sides, you get a suspiciously good fit, but not as suspicious as with 6 sides. But with 8 or 20-sided dice you get a very bad fit, so bad that its probability underflows to 0. This is because the corresponding moduli, 9 and 21, are composite, which means some of the cells in our chi-square test will have no observations. (Suppose m has a proper factor a. Then if a prime p were congruent to a mod m, p would be have to be divisible by a.)

Update: See the next post for a systematic look at different moduli.

You don’t have to use “dice” that correspond to regular solids. You could consider 10-sided “dice,” for example. For such numbers it may be easier to think of spinners than dice, a spinner with 10 equal arc segments it could fall into.

Related post: Probability that a number is prime

Click to learn more about Bayesian statistics consulting