Order Statistics and Benford's Law

Abstract

Fix a base 𝐵>1 and let 𝜁 have the standard exponential distribution; the distribution of digits of 𝜁 base 𝐵 is known to be very close to Benford's law. If there exists a 𝐶 such that the distribution of digits of 𝐶 times the elements of some set is the same as that of 𝜁, we say that set exhibits shifted exponential behavior base 𝐵. Let 𝑋₁,…,𝑋_𝑁 be i.i.d.r.v. If the 𝑋_𝑖's are Unif, then as 𝑁→∞ the distribution of the digits of the differences between adjacent order statistics converges to shifted exponential behavior. If instead 𝑋_𝑖's come from a compactly supported distribution with uniformly bounded first and second derivatives and a second-order Taylor series expansion at each point, then the distribution of digits of any 𝑁^𝛿 consecutive differences and all 𝑁−1 normalized differences of the order statistics exhibit shifted exponential behavior. We derive conditions on the probability density which determine whether or not the distribution of the digits of all the unnormalized differences converges to Benford's law, shifted exponential behavior, or oscillates between the two, and show that the Pareto distribution leads to oscillating behavior.

1. Introduction

Benford's law gives the expected frequencies of the digits in many tabulated data. It was first observed by Newcomb in the 1880s, who noticed that pages of numbers starting with a 1 in logarithm tables were significantly more worn than those starting with a 9. In 1938 Benford [1] observed the same digit bias in a variety of phenomena. From his observations, he postulated that in many datasets, more numbers began with a 1 than with a 9; his investigations (with 20,229 observations) supported his belief. See [2, 3] for a description and history, and [4] for an extensive bibliography.

For any base 𝐵>1, we may uniquely write a positive 𝑥∈ℝ as 𝑥=𝑀_𝐵(𝑥)⋅𝐵^𝑘, where 𝑘∈ℤ and 𝑀_𝐵(𝑥) (called the mantissa) is in [1,𝐵). A sequence of positive numbers {𝑎_𝑛} is Benford base𝐵 if the probability of observing a mantissa of 𝑎_𝑛 base 𝐵 of at most 𝑠 is log_𝐵𝑠. More precisely, for 𝑠∈[1,𝐵], we havelim_𝑁→∞#{𝑛≤𝑁∶1≤𝑀_𝐵(𝛼_𝑛)≤𝑠}𝑁=log_𝐵𝑠. (1.1) Benford behavior for continuous functions is defined analogously. (If the functions are not positive, we study the distribution of the digits of the absolute value of the function.) Thus, working base 10 we find the probability of observing a first the probability of observing a first digit of 𝑑 is log₁₀(𝑑+1)−log₁₀(𝑑), implying that about 30% of the time the first digit is a 1.

We can prove many mathematical systems follow Benford's law, ranging from recurrence relations [5] to 𝑛! [6] to iterates of power, exponential and rational maps, as well as Newton's method [7–9]; to chains of random variables and hierarchical Bayesian models [10]; to values of 𝐿-functions near the critical line; to characteristic polynomials of random matrix ensembles and iterates of the 3𝑥+1-Map [11, 12]; as well as to products of random variables [13]. We also see Benford's law in a variety of natural systems, such as atomic physics [14], biology [15], and geology [16]. Applications of Benford's law range from rounding errors in computer calculations (see [17, page 255]) to detecting tax (see [18, 19]) and voter fraud (see [20]).

This work is motivated by two observations (see Remark 1.9 for more details). First, since Benford's seminal paper, many investigations have shown that amalgamating data from different sources leads to Benford behavior; second, many standard probability distributions are close to Benford behavior. We investigate the distribution of digits of differences of adjacent ordered random variables. For any 𝛿<1, if we study at most 𝑁^𝛿 consecutive differences of a dataset of size 𝑁, the resulting distribution of leading digits depends very weakly on the underlying distribution of the data, and closely approximates Benford's law. We then investigate whether or not studying all the differences leads to Benford behavior; this question is inspired by the first observation above, and has led to new tests for data integrity (see [21]). These tests are quick and easy-to-apply, and have successfully detected problems with some datasets, thus providing a practical application of our main results.

Proving our results requires analyzing the distribution of digits of independent random variables drawn from the standard exponential, and quantifying how close the distribution of digits of a random variable with the standard exponential distribution is to Benford's law. Leemis et al. [22] have observed that the standard exponential is quite close to Benford's law; this was proved by Engel and Leuenberger [23], who showed that the maximum difference in the cumulative distribution function from Benford's law (base 10) is at least .029 and at most .03. We provide an alternate proof of this result in the appendix using a different technique, as well as showing that there is no base 𝐵 such that the standard exponential distribution is Benford base 𝐵 (Corollary A.2).

Both proofs apply Fourier analysis to periodic functions. In [23, equation (5)], the main step is interchanging an integration and a limit. Our proof is based on applying Poisson summation to the derivative of the cumulative distribution function of the logarithms modulo 1, 𝐹_𝐵. Benford's law is equivalent to 𝐹_𝐵(𝑏)=𝑏, which by calculus is the same as 𝐹^′_𝐵(𝑏)=1 and 𝐹_𝐵(0)=0. Thus, studying the deviation of 𝐹^′_𝐵(𝑏) from 1 is a natural way to investigate the deviations from Benford behavior. We hope the details of these calculations may be of use to others in investigating related problems (Poisson summation has been fruitfully used by Kontorovich and Miller [11] and Jang et al. [10] in proving many systems are Benford's; see also [24]).

1.1. Definitions

A sequence {𝑎_𝑛}^∞_𝑛=1⊂[0,1] is equidistributed iflim_𝑁→∞#{𝑛∶𝑛≤𝑁,𝑎_𝑛∈[𝑎,𝑏]}𝑁=𝑏−𝑎∀[𝑎,𝑏]⊂[0,1]. (1.2)Similarly, a continuous random variable on [0,∞), whose probability density function is 𝑝, is equidistributed modulo 1 iflim_𝑇→∞∫^𝑇₀𝜒_𝑎,𝑏(𝑥)𝑝(𝑥)𝑑𝑥∫^𝑇₀𝑝(𝑥)𝑑𝑥=𝑏−𝑎, (1.3)for any [𝑎,𝑏]⊂[0,1], where 𝜒_𝑎,𝑏(𝑥)=1 for 𝑥mod1∈[𝑎,𝑏] and 0 otherwise.

A positive sequence (or values of a function) is Benford base 𝐵 if and only if its base 𝐵 logarithms are equidistributed modulo 1; this equivalence is at the heart of many investigations of Benford's law (see [6, 25] for a proof).

We use the following notations for the various error terms.

(1) Let ℰ(𝑥) denote an error of at most 𝑥 in absolute value; thus 𝑓(𝑏)=𝑔(𝑏)+ℰ(𝑥) means |𝑓(𝑏)−𝑔(𝑏)|≤𝑥.(2) Big-Oh notation: for 𝑔(𝑥) a nonnegative function, we say 𝑓(𝑥)=𝑂(𝑔(𝑥)) if there exist an 𝑥₀ and a 𝐶>0 such that, for all 𝑥≥𝑥₀, |𝑓(𝑥)|≤𝐶𝑔(𝑥).

The following theorem is the starting point for investigating the distribution of digits of order statistics.

Theorem 1.1. Let 𝜁 have the standard (unit) exponential distributionProb(𝜁∈[𝛼,𝛽])=^𝛽_𝛼𝑒^−𝑡𝑑𝑡,[𝛼,𝛽]∈[0,∞). (1.4)For 𝑏∈[0,1], let 𝐹_𝐵(𝑏) be the cumulative distribution function of log_𝐵𝜁mod1; thus 𝐹_𝐵(𝑏)∶=Prob(log_𝐵𝜁mod1∈[0,𝑏]). Then, for all 𝑀≥2,𝐹^′_𝐵(𝑏)=1+2^∞_𝑚=1Re(𝑒^{−2𝜋𝑖𝑚𝑏}Γ(1+2𝜋𝑖𝑚log𝐵))=1+2^𝑀−1_𝑚=1Re(𝑒^{−2𝜋𝑖𝑚𝑏}Γ(1+2𝜋𝑖𝑚√log𝐵))+ℰ(42𝜋𝑐₁(𝐵)𝑒^{−(𝜋2−𝑐2(𝐵))𝑀/log𝐵}), (1.5)where 𝑐₁(𝐵),𝑐₂(𝐵) are constants such that for all 𝑚≥𝑀≥2, one has 𝑒^{2𝜋2𝑚/log𝐵}−𝑒^{−2𝜋2𝑚/log𝐵}≥𝑒^{2𝜋2𝑚/log𝐵}𝑐²₁,𝑚(𝐵)log𝐵≤𝑒^{2𝑐2(𝐵)𝑚/log𝐵},1−𝑒^{(𝜋2−𝑐2(𝐵))𝑀/log𝐵}≥1√2. (1.6)For 𝐵∈[𝑒,10], take 𝑐₁(𝐵)=√2 and 𝑐₂(𝐵)=1/5, which giveProb(log𝜁mod1∈[𝑎,𝑏])=𝑏−𝑎+2𝑟𝜋⋅sin(𝜋(𝑏+𝑎)+𝜃)⋅sin(𝜋(𝑏−𝑎))+ℰ(6.32⋅10⁻⁷), (1.7)with 𝑟≈0.000324986, 𝜃≈1.32427186, andProb(log₁₀𝜁mod1∈[𝑎,𝑏])=𝑏−𝑎+2𝑟₁𝜋sin(𝜋(𝑏+𝑎)−𝜃₁−𝑟)⋅sin(𝜋(𝑏−𝑎))₂𝜋sin(2𝜋(𝑏+𝑎)+𝜃₂)⋅sin(2𝜋(𝑏−𝑎))+ℰ(8.5⋅10⁻⁵) (1.8)with𝑟₁≈0.0569573,𝜃₁𝑟≈0.8055888,₂≈0.0011080,𝜃₂≈0.1384410. (1.9)

The above theorem was proved in [23]; we provide an alternate proof in Appendix A. As remarked earlier, our technique consists of applying Poisson summation to the derivative of the cumulative distribution function of the logarithms modulo 1; it is then very natural and easy to compare deviations from the resulting distribution and the uniform distribution (if a dataset satisfies Benford's law, then the distribution of its logarithms is uniform). Our series expansions are obtained by applying properties of the Gamma function.

Definition 1.2 (Definition exponential behavior, shifted exponential behavior). Let 𝜁 have the standard exponential distribution, and fix a base 𝐵. If the distribution of the digits of a set is the same as the distribution of the digits of 𝜁, then one says that the set exhibits exponential behavior (base 𝐵). If there is a constant 𝐶>0 such that the distribution of digits of all elements multiplied by 𝐶 is exponential behavior, then one says that the system exhibits shifted exponential behavior (with shift of log_𝐵𝐶mod1).

We briefly describe the reasons behind this notation. One important property of Benford's law is that it is invariant under rescaling; many authors have used this property to characterize Benford behavior. Thus, if a dataset is Benford base 𝐵, and we fix a positive number 𝐶, so is the dataset obtained by multiplying each element by 𝐶. This is clear if, instead of looking at the distribution of the digits, we study the distribution of the base 𝐵 logarithms modulo 1. Benford's law is equivalent to the logarithms modulo 1 being uniformly distributed (see, e.g., [6, 25]); the effect of multiplying all entries by a fixed constant simply translates the uniform distribution modulo 1, which is again the uniform distribution.

The situation is different for exponential behavior. Multiplying all elements by a fixed constant 𝐶 (where 𝐶≠𝐵^𝑘 for some 𝑘∈ℤ) does not preserve exponential behavior; however, the effect is easy-to-describe. Again looking at the logarithms, exponential behavior is equivalent to the base 𝐵 logarithms modulo 1 having a specific distribution which is almost equal to the uniform distribution (at least if the base 𝐵 is not too large). Multiplying by a fixed constant 𝐶≠𝐵^𝑘 shifts the logarithm distribution by log_𝐵𝐶mod1.

1.2. Results for Differences of Orders Statistics

We consider a simple case first, and show how the more general case follows. Let 𝑋₁,…,𝑋_𝑁 be independent identically distributed from the uniform distribution on [0,𝐿]. We consider 𝐿 fixed and study the limit as 𝑁→∞. Let 𝑋_1∶𝑁,…,𝑋_𝑁∶𝑁 be the 𝑋_𝑖's in increasing order. The 𝑋_𝑖∶𝑁 are called the order statistics, and satisfy 0≤𝑋_1∶𝑁≤𝑋_2∶𝑁≤⋯≤𝑋_𝑁∶𝑁≤𝐿. We investigate the distribution of the leading digits of the differences between adjacent 𝑋_𝑖∶𝑁's, 𝑋_{𝑖+1∶𝑁}−𝑋_𝑖∶𝑁. For convenience, we periodically continue the data and set 𝑋_{𝑖+𝑁∶𝑁}=𝑋_𝑖∶𝑁+𝐿. As we have 𝑁 differences in an interval of size 𝐿, the average value of 𝑋_{𝑖+1∶𝑁}−𝑋_𝑖∶𝑁 is of size 𝐿/𝑁, and it is sometimes easier to study the normalized differences𝑍_𝑖;𝑁=𝑋_{𝑖+1∶𝑁}−𝑋_𝑖∶𝑁𝐿/𝑁. (1.10)As the 𝑋_𝑖's are drawn from a uniform distribution, it is a standard result that as 𝑁→∞, the 𝑍_𝑖;𝑁's are independent random variables, each having the standard exponential distribution. Thus, as 𝑁→∞, the probability that 𝑍_𝑖;𝑁∈[𝑎,𝑏] tends to ∫^𝑏_𝑎𝑒^−𝑡𝑑𝑡 (see [26, 27] for proofs).

For uniformly distributed random variables, if we know the distribution of log_𝐵𝑍_𝑖;𝑁mod1, then we can immediately determine the distribution of the digits of the 𝑋_{𝑖+1∶𝑁}−𝑋_𝑖∶𝑁 base 𝐵 becauselog_𝐵𝑍_𝑖;𝑁=log_𝐵(𝑋_{𝑖+1∶𝑁}−𝑋_𝑖∶𝑁𝐿/𝑁)=log_𝐵(𝑋_{𝑖+1∶𝑁}−𝑋_𝑖∶𝑁)−log_𝐵(𝐿𝑁). (1.11)

As 𝑍_𝑖;𝑁 are independent with the standard exponential distribution as 𝑁→∞; if 𝑋_𝑖 are independent uniformly distributed, the behavior of the digits of the differences 𝑋_{𝑖+1∶𝑁}−𝑋_𝑖∶𝑁 is an immediate consequence of Theorem 1.1. Theorem 1.3 (Shifted exponential behavior of differences of independent uniformly distributed random variables). Let 𝑋₁,…,𝑋_𝑁 be independently distributed from the uniform distribution on [0,𝐿], and let 𝑋_1∶𝑁,…,𝑋_𝑁∶𝑁 be 𝑋_𝑖's in an increasing order. As 𝑁→∞, the distribution of the digits (base 𝐵) of the differences 𝑋_{𝑖+1∶𝑁}−𝑋_𝑖∶𝑁 converges to shifted exponential behavior, with a shift of log_𝐵(𝐿/𝑁)mod1.

A similar result holds for other distributions. Theorem 1.4 (Shifted exponential behavior of subsets of differences of independent random variables). Let 𝑋₁,…,𝑋_𝑁 be independent, identically distributed random variables whose density 𝑓(𝑥) has a second-order Taylor series at each point with first and second derivatives uniformly bounded, and let the 𝑋_𝑖∶𝑁's be the 𝑋_𝑖's in increasing order. Fix a 𝛿∈(0,1). Then, as 𝑁→∞ the distribution of the digits (base 𝐵) of 𝑁^𝛿 consecutive differences 𝑋_{𝑖+1∶𝑁}−𝑋_𝑖∶𝑁 converges to shifted exponential behavior, provided that 𝑋_𝑖∶𝑁's are from a region where 𝑓(𝑥) is nonzero.

The key ingredient in this generalization is that the techniques, which show that the differences between uniformly distributed random variables become independent exponentially distributed random variables, can be modified to handle more general distributions.

We restricted ourselves to a subset of all consecutive spacings because the normalization factor changes throughout the domain. The shift in the shifted exponential behavior depends on which set of 𝑁^𝛿 differences we study, coming from the variations in the normalizing factors. Within a bin of 𝑁^𝛿 differences, the normalization factor is basically constant, and we may approximate our density with a uniform distribution. It is possible for these variations to cancel and yield Benford behavior for the digits of all the unnormalized differences. Such a result is consistent with the belief that amalgamation of data from many different distributions becomes Benford; however, this is not always the case (see Remark 1.6). From Theorems 1.1 and 1.4, we obtain the following theorem. Theorem 1.5 (Benford behavior for all the differences of independent random variables). Let 𝑋₁,…,𝑋_𝑁 be independent identically distributed random variables whose density 𝑓(𝑥) is compactly supported and has a second-order Taylor series at each point with first and second derivatives uniformly bounded. Let the 𝑋_𝑖∶𝑁's be the 𝑋_𝑖's in an increasing order 𝐹(𝑥) be the cumulative distribution function for 𝑓(𝑥), and fix a 𝛿∈(0,1). Let 𝐼(𝜖,𝛿,𝑁)=[𝜖𝑁^1−𝛿,𝑁^1−𝛿−𝜖𝑁^1−𝛿]. For each fixed 𝜖∈(0,1/2), assume that
(i)𝑓(𝐹⁻¹(𝑘𝑁^𝛿−1) is not too small for 𝑘∈𝐼(𝜖,𝛿,𝑁)lim_𝑁→∞max_{𝑘∈𝐼(𝜖,𝛿,𝑁)}min(𝑁^{−(𝜖+𝛿/2)},𝑁^𝛿−1)𝑓(𝐹⁻¹(𝑘𝑁^𝛿−1))=0, (1.12)(ii)log_𝐵𝑓(𝐹⁻¹(𝑘𝑁^𝛿−1)mod1 is equidistributed: for all [𝛼,𝛽]⊂[0,1]lim_𝑁→∞#{𝑘∈𝐼(𝜖,𝛿,𝑁)∶log_𝐵𝑓(𝐹⁻¹(𝑘𝑁^𝛿−1))mod1∈[𝛼,𝛽]}𝑁^𝛿=𝛽−𝛼. (1.13) Then, if 𝜖>max(0,1/3−𝛿/2) and 𝜖<𝛿/2, the distribution of the digits of the 𝑁−1 differences 𝑋_{𝑖+1∶𝑁}−𝑋_𝑖∶𝑁 converges to Benford's law (base 𝐵) as 𝑁→∞.

Remark 1.6. The conditions of Theorem 1.5 are usually not satisfied. We are unaware of any situation where (1.13) holds; we have included Theorem 1.5 to give a sufficient condition of what is required to have Benford's law satisfied exactly, and not just approximately. In Lemma 3.3, we show with: Example 3.3 shows that the conditions fail for the Pareto distribution, and the limiting behavior oscillates between Benford and a sum of shifted exponential behavior. (If several datasets each exhibit shifted exponential behavior but with distinct shifts, then the amalgamated dataset is closer to Benford's law than any of the original datasets. This is apparent by studying the logarithms modulo 1. The differences between these densities and Benford's law will look like Figure 1(b) (except, of course, that different shifts will result in shifting the plot modulo 1). The key observation is that the unequal shifts mean that we do not have reinforcements from the peaks of modulo 1 densities being aligned, and thus the amalgamation will decrease the maximum deviations.) The arguments generalize to many densities whose cumulative distribution functions have tractable closed-form expressions (e.g., exponential, Weibull, or 𝑓(𝑥)=𝑒^−𝑒𝑥𝑒^𝑥).

(a)

(b)

Figure 1 

All 499 999 differences of adjacent order statistics from 500 000 independent random variables from the Pareto distribution with minimum value and variance 1. (a) Observed digits of scaled differences of adjacent random variables versus Benford's law; (b) scaled observed minus Benford's law (cumulative distribution of base 10 logarithms).

The situation is very different if instead we study normalized differences:^0𝑥0𝑍^02𝑑𝑐_𝑖∶𝑁=𝑋_{𝑖+1∶𝑁}−𝑋_𝑖∶𝑁1/𝑁𝑓(𝑋_𝑖∶𝑁), (1.14)note if 𝑓(𝑥)=1/𝐿 is the uniform distribution on [0,𝐿], (1.14) reduces to (1.10). Theorem 1.7 (Shifted exponential behavior for all the normalized differences of independent random variables). Assume the probability distribution 𝑓 satisfies the conditions of Theorem 1.5 and (1.12) and ^0𝑥0𝑍^02𝑑𝑐_𝑖;𝑁 is as in (1.14). Then, as 𝑁→∞, the distribution of the digits of ^0𝑥0𝑍^02𝑑𝑐_𝑖∶𝑁 converges to shifted exponential behavior.

Remark 1.8. Appropriately scaled, the distribution of the digits of the differences is universal, and is the exponential behavior of Theorem 1.1. Thus, Theorem 1.7 implies that the natural quantity to study is the normalized differences of the order statistics, not the differences (see also Remark 3.5). With additional work, we could study densities with unbounded support and show that, through truncation, we can get arbitrarily close to shifted exponential behavior.

Remark 1.9. The main motivation for this work is the need for improved ways of assessing the authenticity and integrity of scientific and corporate data. Benford's law has been successfully applied to detecting income tax, corporate, and voter fraud (see [18–20]); in [21], we use these results to derive new statistical tests to examine data authenticity and integrity. Early applications of these tests to financial data showed that it could detect errors in data downloads, rounded data, and inaccurate ordering of data. These attributes are not easily observable from an analysis of descriptive statistics, and detecting these errors can help managers avoid costly decisions based on erroneous data.

The paper is organized as follows. We prove Theorem 1.1 in Appendix A by using Poisson summation to analyze 𝐹^′_𝐵(𝑏). Theorem 1.3 follows from the results of the order statistics of independent uniform variables. The proof of Theorem 1.4 is similar, and given in Section 2. In Section 3, we prove Theorems 1.5 and 1.7.

2. Proofs of Theorems 1.3 and 1.4

Theorem 1.3 is a consequence of the fact that the normalized differences between the order statistics drawn from the uniform distribution converge to being independent standard exponentials. The proof of Theorem 1.4 proceeds similarly. Specifically, over a short enough region, any distribution with a second-order Taylor series at each point with first and second derivatives uniformly bounded is well approximated by a uniform distribution.

To prove Theorem 1.4, it suffices to show that if 𝑋₁,…,𝑋_𝑁 are drawn from a sufficiently nice distribution, then for any fixed 𝛿∈(0,1) the limiting behavior of the order statistics of 𝑁^𝛿 adjacent 𝑋_𝑖's becomes Poissonian (i.e., the 𝑁^𝛿−1 normalized differences converge to being independently distributed from the standard exponential). We prove this below for compactly supported distributions 𝑓(𝑥) that have a second-order Taylor series at each point with the first and second derivatives uniformly bounded, and when the 𝑁^𝛿 adjacent 𝑋_𝑖's are from a region where 𝑓(𝑥) is bounded away from zero.

For each 𝑁, consider intervals [𝑎_𝑁,𝑏_𝑁] such that ∫^𝑏𝑁_𝑎𝑁𝑓(𝑥)𝑑𝑥=𝑁^𝛿/𝑁; thus, the proportion of the total mass in such intervals is 𝑁^𝛿−1. We fix such an interval for our arguments. For each 𝑖∈{1,…,𝑁}, let𝑤_𝑖=1,if𝑋_𝑖∈[𝑎_𝑁,𝑏_𝑁]0,otherwise. (2.1) Note 𝑤_𝑖 is 1 with probability 𝑁^𝛿−1 and 0 with probability 1−𝑁^𝛿−1; 𝑤_𝑖 is a binary indicator random variable, telling us whether or not 𝑋_𝑖∈[𝑎_𝑁,𝑏_𝑁]. Thus,𝔼[^𝑁_𝑖=1𝑤_𝑖]=𝑁^𝛿,Var(^𝑁_𝑖=1𝑤_𝑖)=𝑁^𝛿⋅(1−𝑁^𝛿−1). (2.2)Let 𝑀_𝑁 be the number of 𝑋_𝑖 in [𝑎_𝑁,𝑏_𝑁], and let 𝛽_𝑁 be any nondecreasing sequence tending to infinity (in the course of the proof, we will find that we may take any sequence with 𝛽_𝑁=𝑂(𝑁^𝛿/2)). By (2.2) and the central limit theorem (which we may use as 𝑤_𝑖's satisfy the Lyapunov condition), with probability tending to 1, we have𝑀_𝑁=𝑁^𝛿+𝑂(𝛽_𝑁𝑁^𝛿/2). (2.3)

We assume that in the interval [𝑎_𝑁,𝑏_𝑁] there exist constants 𝑐 and 𝐶 such that whenever 𝑥∈[𝑎_𝑁,𝑏_𝑁], 0<𝑐<𝑓(𝑥)<𝐶<∞; we assume that these constants hold for all regions investigated and for all 𝑁. (If our distribution has unbounded support, for any 𝜖>0, we can truncate it on both sides so that the omitted probability is at most 𝜖. Our result is then trivially modified to be within 𝜖 of shifted exponential behavior.) Thus,𝑐⋅(𝑏_𝑁−𝑎_𝑁)≤^𝑏𝑁_𝑎𝑁𝑓(𝑥)𝑑𝑥=𝑁^𝛿−1≤𝐶(𝑏_𝑁−𝑎_𝑁), (2.4)implying that 𝑏_𝑁−𝑎_𝑁 is of size 𝑁^𝛿−1. If we assume that 𝑓(𝑥) has at least a second-order Taylor expansion, then𝑓(𝑥)=𝑓(𝑎_𝑁)+𝑓^′(𝑎_𝑁)(𝑥−𝑎_𝑁)+𝑂((𝑥−𝑎_𝑁)²)=𝑓(𝑎_𝑁)+𝑓^′(𝑎_𝑁)(𝑥−𝑎_𝑁)+𝑂(𝑁^2𝛿−2). (2.5) As we assume that the first and second derivatives are uniformly bounded, as well as 𝑓 being bounded away from zero in the intervals under consideration, all Big-Oh constants below are independent of 𝑁. Thus,𝑏_𝑁−𝑎_𝑁=𝑁^𝛿−1𝑓(𝑎_𝑁)+𝑂(𝑁^2𝛿−2). (2.6)

We now investigate the order statistics of the 𝑀_𝑁 of the 𝑋_𝑖's that lie in [𝑎_𝑁,𝑏_𝑁]. We know ∫^𝑏𝑁_𝑎𝑁𝑓(𝑥)𝑑𝑥=𝑁^𝛿−1; by setting 𝑔_𝑁(𝑥)=𝑓(𝑥)𝑁^1−𝛿, then 𝑔_𝑁(𝑥) is the conditional density function for 𝑋_𝑖, given that 𝑋_𝑖∈[𝑎_𝑁,𝑏_𝑁]. Thus, 𝑔_𝑁(𝑥) integrates to 1, and for 𝑥∈[𝑎_𝑁,𝑏_𝑁], we have𝑔_𝑁(𝑥)=𝑓(𝑎_𝑁)⋅𝑁^1−𝛿+𝑓^′(𝑎_𝑁)(𝑥−𝑎_𝑁)⋅𝑁^1−𝛿+𝑂(𝑁^𝛿−1). (2.7)

We have an interval of size 𝑁^𝛿−1/𝑓(𝑎_𝑁)+𝑂(𝑁^2𝛿−2), and 𝑀_𝑁=𝑁^𝛿+𝑂(𝛽_𝑁𝑁^𝛿/2) of the 𝑋_𝑖 lying in the interval (remember that 𝛽_𝑁 are any nondecreasing sequence tending to infinity). Thus, with probability tending to 1, the average spacing between adjacent ordered 𝑋_𝑖 is𝑁^𝛿−1/𝑓(𝑎_𝑁)+𝑂(𝑁^2𝛿−2)𝑀_𝑁=(𝑓(𝑎_𝑁)𝑁)⁻¹+𝑁⁻¹⋅𝑂(𝛽_𝑁𝑁^−𝛿/2+𝑁^𝛿−1), (2.8)in particular, we see that we must choose 𝛽_𝑁=𝑂(𝑁^𝛿/2). As 𝛿∈(0,1), if we fix a 𝑘 such that 𝑋_𝑘∈[𝑎_𝑁,𝑏_𝑁], then we expect the next 𝑋_𝑖 to the right of 𝑋_𝑘 to be about 𝑡/𝑁𝑓(𝑎_𝑁) units away, where 𝑡 is of size 1. For a given 𝑋_𝑘, we can compute the conditional probability that the next 𝑋_𝑖 is between 𝑡/𝑁𝑓(𝑎_𝑁) and (𝑡+Δ𝑡)/𝑁𝑓(𝑎_𝑁) units to the right. It is simply the difference of the probability that all the other 𝑀_𝑁−1 of the 𝑋_𝑖's in [𝑎_𝑁,𝑏_𝑁] are not in the interval [𝑋_𝑘,𝑋_𝑘+𝑡/𝑁𝑓(𝑎_𝑁)] and the probability that all other 𝑋_𝑖 in [𝑎_𝑁,𝑏_𝑁] are not in the interval [𝑋_𝑘,𝑋_𝑘+(𝑡+Δ𝑡)/𝑁𝑓(𝑎_𝑁)]; note that we are using the wrapped interval [𝑎_𝑁,𝑏_𝑁].

Some care is required in these calculations. We have a conditional probability as we assume that 𝑋_𝑘∈[𝑎_𝑁,𝑏_𝑁] and that exactly 𝑀_𝑁 of the 𝑋_𝑖 are in [𝑎_𝑁,𝑏_𝑁]. Thus, these probabilities depend on two random variables, namely, 𝑋_𝑘 and 𝑀_𝑁. This is not a problem in practice, however (e.g., 𝑀_𝑁 is tightly concentrated about its mean value).

Recalling our expansion for 𝑔_𝑁(𝑥) (and that 𝑏_𝑁−𝑎_𝑁=𝑁^𝛿−1/𝑓(𝑎_𝑁)+𝑂(𝑁^2𝛿−2) and 𝑡 is of size 1), after simple algebra, we find that with probability tending to 1, for a given 𝑋_𝑘 and 𝑀_𝑁, the first probability is(1−^{𝑋𝑘+𝑡/𝑁𝑓(𝑎𝑁)}_𝑋𝑘𝑔_𝑁(𝑥)𝑑𝑥)^𝑀𝑁−1. (2.9)The above integral equals 𝑡𝑁^𝛿+𝑂(𝑁⁻¹) (use the Taylor series expansion in (2.7) and note that the interval [𝑎_𝑁,𝑏_𝑁] is of size 𝑂(𝑁^𝛿−1)). Using (2.3), it is easy to see that this is a.s. equal to(1−𝑡+𝑂(𝑁^𝛿−1+𝛽_𝑁𝑁^−𝛿/2)𝑀_𝑁)^𝑀𝑁−1. (2.10)We, therefore, find that as 𝑁→∞, the probability that 𝑀_𝑁−1 of the 𝑋_𝑖's (𝑖≠𝑘) are in [𝑎_𝑁,𝑏_𝑁]⧵[𝑋_𝑘,𝑋_𝑘+𝑡/𝑁𝑓(𝑎_𝑁)], conditioned on 𝑋_𝑘 and 𝑀_𝑁, converges to 𝑒^−𝑡. (Some care is required, as the exceptional set in our a.s. statement can depend on 𝑡. This can be surmounted by taking expectations with respect to our conditional probabilities and applying the dominated convergence theorem.)

The calculation of the second probability, the conditional probability that the 𝑀_𝑁−1 other 𝑋_𝑖' that are [𝑎_𝑁,𝑏_𝑁] not in the interval [𝑋_𝑘,𝑋_𝑘+(𝑡+Δ𝑡)/𝑁𝑓(𝑎_𝑁)], given 𝑋_𝑘 and 𝑀_𝑁, follows analogously by replacing 𝑡 with 𝑡+Δ𝑡 in the previous argument. We thus find that this probability is 𝑒^{−(𝑡+Δ𝑡)}. As_𝑡^𝑡+Δ𝑡𝑒^−𝑢𝑑𝑢=𝑒^−𝑡−𝑒^{−(𝑡+Δ𝑡)}, (2.11)we find that the density of the difference between adjacent order statistics tends to the standard (unit) exponential density; thus, the proof of Theorem 1.4 now follows from Theorem 1.3.

3. Proofs of Theorems 1.5 and 1.7

We generalize the notation from Section 2. Let 𝑓(𝑥) be any distribution with a second-order Taylor series at each point with first and second derivatives uniformly bounded, and let 𝑋_1∶𝑁,…,𝑋_𝑁∶𝑁 be the order statistics. We fix a 𝛿∈(0,1), and for 𝑘∈{1,…,𝑁^1−𝛿}, we consider bins [𝑎_𝑘;𝑁,𝑏_𝑘;𝑁] such that^{𝑏𝑘;𝑁}_{𝑎𝑘;𝑁}𝑁𝑓(𝑥)𝑑𝑥=^𝛿𝑁=𝑁^𝛿−1, (3.1)there are 𝑁^1−𝛿 such bins. By the central limit theorem (see (2.3)), if 𝑀_𝑘;𝑁 is the number of order statistics in [𝑎_𝑘;𝑁,𝑏_𝑘;𝑁], then, provided that 𝜖>max(0,1/3−𝛿/2) with probability tending to 1, we have𝑀_𝑘;𝑁=𝑁^𝛿+𝑂(𝑁^𝜖+𝛿/2), (3.2)of course we also require 𝜖<𝛿/2, as, otherwise, the error term is larger than the main term.

Remark 3.1. Before we considered just one fixed interval; as we are studying 𝑁^1−𝛿 intervals simultaneously, we need 𝜖 in the exponent so that with high probability, all intervals have to first order 𝑁^𝛿 order statistics. For the arguments below, it would have sufficed to have an error of size 𝑂(𝑁^𝛿−𝜖). We thank the referee for pointing out that 𝜖>1/3−𝛿/2, and provide his argument in Appendix B.

Similar to (2.8), the average spacing between adjacent order statistics in [𝑎_𝑘;𝑁,𝑏_𝑘;𝑁] is(𝑓(𝑎_𝑘;𝑁)𝑁)⁻¹+𝑁⁻¹⋅𝑂(𝑁^{−(𝜖+𝛿/2)}+𝑁^𝛿−1). (3.3)Note that (3.3) is the generalization of (1.11); if 𝑓 is the uniform distribution on [0,𝐿], then 𝑓(𝑎_𝑘;𝑁)=1/𝐿. By Theorem 1.4, as 𝑁→∞, the distribution of digits of the differences in each bin converges to shifted exponential behavior; however, the variation in the average spacing between bins leads to bin-dependent shifts in the shifted exponential behavior.

Similar to (1.11), we can study the distribution of digits of the differences of the normalized order statistics. If 𝑋_𝑖∶𝑁 and 𝑋_{𝑖+1∶𝑁} are in [𝑎_𝑘;𝑁,𝑏_𝑘;𝑁], then𝑍_𝑖;𝑁=𝑋_{𝑖+1∶𝑁}−𝑋_𝑖∶𝑁(𝑓(𝑎_𝑘;𝑁)𝑁)⁻¹+𝑁⁻¹⋅𝑂(𝑁^{−(𝜖+𝛿/2)}+𝑁^𝛿−1),log_𝐵𝑍_𝑖;𝑁=log_𝐵(𝑋_{𝑖+1∶𝑁}−𝑋_𝑖∶𝑁)+log_𝐵𝑁−log_𝐵(𝑓(𝑎_𝑘;𝑁)⁻¹+𝑂(𝑁^{−(𝜖+𝛿/2)}+𝑁^𝛿−1)). (3.4)Note we are using the same normalization factor for all differences between adjacent order statistics in a bin. Later, we show that we may replace 𝑓(𝑎_𝑘;𝑁) with 𝑓(𝑋_𝑖∶𝑁). As we study all 𝑋_{𝑖+1∶𝑁}−𝑋_𝑖∶𝑁 in the bin [𝑎_𝑘;𝑁,𝑏_𝑘;𝑁], it is useful to rewrite the above aslog_𝐵(𝑋_{𝑖+1∶𝑁}−𝑋_𝑖∶𝑁)=log_𝐵𝑍_𝑖;𝑁−log_𝐵𝑁+log_𝐵(𝑓(𝑎_𝑘;𝑁)⁻¹+𝑂(𝑁^{−(𝜖+𝛿/2)}+𝑁^𝛿−1)). (3.5)We have 𝑁^1−𝛿 bins, so 𝑘∈{1,…,𝑁^1−𝛿}. As we only care about the limiting behavior, we may safely ignore the first and last bins. We may, therefore, assume that each 𝑎_𝑘;𝑁 is finite, and 𝑎_𝑘+1;𝑁=𝑏_𝑘;𝑁. (Of course, we know that both quantities are finite as we assumed that our distribution has compact support. We remove the last bins to simplify generalizations to noncompactly supported distributions.)

Let 𝐹(𝑥) be the cumulative distribution function for 𝑓(𝑥). Then,𝐹(𝑎_𝑘;𝑁)=(𝑘−1)𝑁^𝛿𝑁=(𝑘−1)𝑁^𝛿−1. (3.6)For notational convenience, we relabel the bins so that 𝑘∈{0,…,𝑁^1−𝛿−1}; thus 𝐹(𝑎_𝑘;𝑁)=𝑘𝑁^𝛿−1.

We now prove our theorems which determine when these bin-dependent shifts cancel (yielding Benford behavior), or reinforce (yielding sums of shifted exponential behavior). Proof of Theorem 1.5. There are approximately 𝑁^𝛿 differences in each bin [𝑎_𝑘;𝑁,𝑏_𝑘;𝑁]. By Theorem 1.4, the distribution of the digits of the differences in each bin converges to shifted exponential behavior. As we assume that the first and second derivatives of 𝑓 are uniformly bounded, the Big-Oh constants in Section 2 are independent of the bins. The shift in the shifted exponential behavior in each bin is controlled by the last two terms on the right-hand side of (3.5). The log_𝐵𝑁 shifts the shifted exponential behavior in each bin equally. The bin-dependent shift is controlled by the final termlog_𝐵(𝑓(𝑎_𝑘;𝑁)⁻¹+𝑂(𝑁^{−(𝜖+𝛿/2)}+𝑁^𝛿−1))=−log_𝐵𝑓(𝑎_𝑘;𝑁)+log_𝐵(1+min(𝑁^{−(𝜖+𝛿/2)},𝑁^𝛿−1)𝑓(𝑎_𝑘;𝑁)). (3.7)
Thus, each of the 𝑁^1−𝛿 bins exhibits shifted exponential behavior, with a bin-dependent shift composed of the two terms in (3.7). By (1.12), 𝑓(𝑎_𝑘;𝑁) are not small compared to min(𝑁^{−(𝜖+𝛿/2)},𝑁^𝛿−1), and hence the second term log_𝐵(1+(min(𝑁^{−(𝜖+𝛿/2)},𝑁^𝛿−1)/𝑓(𝑎_𝑘;𝑁))) is negligible. In particular, this factor depends only very weakly on the bin, and tends to zero as 𝑁→∞.
Thus, the bin-dependent shift in the shifted exponential behavior is approximately −log_𝐵𝑓(𝑎_𝑘;𝑁)=−log_𝐵𝑓(𝐹⁻¹(𝑘𝑁^𝛿−1)). If these shifts are equidistributed modulo 1, then the deviations from Benford behavior cancel, and the shifted exponential behavior of each bin becomes Benford behavior for all the differences.

Remark 3.2. Consider the case when the density is a uniform distribution on some interval. Then, all 𝑓(𝐹⁻¹(𝑘𝑁^𝛿−1)) are equal, and each bin has the same shift in its shifted exponential behavior. These shifts, therefore, reinforce each other, and the distribution of all the differences is also shifted exponential behavior, with the same shift. This is observed in numerical experiments (see Theorem 1.3 for an alternate proof).

We analyze the assumptions of Theorem 1.5. The condition from (1.12) is easy-to-check, and is often satisfied. For example, if the probability density is a finite union of monotonic pieces and is zero only finitely often, then (1.12) holds. This is because for 𝑘∈𝐼(𝜖,𝛿,𝑁), 𝐹⁻¹(𝑘𝑁^𝛿−1)∈[𝐹⁻¹(𝜖),𝐹⁻¹(1−𝜖)], and this is, therefore, independent of 𝑁 (if 𝑓 vanishes finitely often, we need to remove small subintervals from 𝐼(𝜖,𝛿,𝑁), but the analysis proceeds similarly). The only difficulty is basically a probability distribution with intervals of zero probability. Thus, (1.12) is a mild assumption.

If we choose any distribution other than a uniform distribution, then 𝑓(𝑥) is not constant; however, (1.13) does not need to hold (i.e., log_𝐵𝑓(𝑎_𝑘;𝑁)mod1 does not need to be equidistributed as 𝑁→∞). For example, consider a Pareto distribution with minimum value 1 and exponent 𝑎>0. The density is𝑓(𝑥)=𝑎𝑥^−𝑎−1if𝑥≥1,0otherwise. (3.8) The Pareto distribution is known to be useful in modelling natural phenomena, and for appropriate choices of exponents, it yields approximately Benford behavior (see [16]).

Example 3.3. If 𝑓 is a Pareto distribution with minimum value 1 and exponent 𝑎>0, then 𝑓 does not satisfy the second condition of Theorem 1.5, (1.13).
To see this, note that the cumulative distribution function of 𝑓 is 𝐹(𝑥)=1−𝑥^−𝑎. As we only care about the limiting behavior, we need only to study 𝑘∈𝐼(𝜖,𝛿,𝑁)=[𝜖𝑁^1−𝛿,𝑁^1−𝛿−𝜖𝑁^1−𝛿]. Therefore, 𝐹(𝑎_𝑘;𝑁)=𝑘𝑁^𝛿−1 implies that𝑎_𝑘;𝑁=(1−𝑘𝑁^𝛿−1)^−1/𝑎,𝑓(𝑎_𝑘;𝑁)=𝑎(1−𝑘𝑁^𝛿−1)^{(𝑎+1)/𝑎}. (3.9)
The condition from (1.12) is satisfied, namely,lim_𝑁→∞max_{𝑘∈𝐼(𝜖,𝛿,𝑁)}min(𝑁^{−(𝜖+𝛿/2)},𝑁^𝛿−1)𝑓(𝑎_𝑘;𝑁)=lim_𝑁→∞max_{𝑘∈𝐼(𝜖,𝛿,𝑁)}min(𝑁^{−(𝜖+𝛿/2)},𝑁^𝛿−1)𝑎(𝑘𝑁^𝛿−1)^{(𝑎+1)/𝑎}=0, (3.10)as 𝑘 is of size 𝑁^1−𝛿.
Let 𝑗=𝑁^1−𝛿−𝑘∈𝐼(𝜖,𝛿,𝑁). Then, the bin-dependent shifts arelog_𝐵𝑓(𝑎_𝑘;𝑁)=𝑎+1𝑎log_𝐵(1−𝑘𝑁^𝛿−1)+log_𝐵𝑎=𝑎+1𝑎log_𝐵(𝑗𝑁^1−𝛿)+log_𝐵𝑎=log_𝐵(𝑗^{(𝑎+1)/𝑎})+log_𝐵(𝑎𝑁^{(1−𝛿)(𝑎+1)/𝑎}). (3.11)Thus, for a Pareto distribution with exponent 𝑎, the distribution of all the differences becomes Benford if and only if 𝑗^{(𝑎+1)/𝑎} is Benford. This follows from the fact that a sequence is Benford if and only if its logarithms are equidistributed. For fixed 𝑚, 𝑗^𝑚 is not Benford (e.g., [6]), and thus the condition from (1.13) fails.

Remark 3.4. We chose to study a Pareto distribution because the distribution of digits of a random variable drawn from a Pareto distribution converges to Benford behavior (base 10) as 𝑎→1; however, the digits of the differences do not tend to Benford (or shifted exponential) behavior. A similar analysis holds for many distributions with good closed-form expressions for the cumulative distribution function. In particular, if 𝑓 is the density of an exponential or Weibull distribution (or 𝑓(𝑥)=𝑒^−𝑒𝑥𝑒^𝑥), then 𝑓 does not satisfy the second condition of Theorem 1.5, (1.13).

Modifying the proof of Theorem 1.5 yields our result on the distribution of digits of the normalized differences.

Proof of Theorem 1.7. If 𝑓 is the uniform distribution, there is nothing to prove. For general 𝑓, rescaling the differences eliminates the bin-dependent shifts. Let^0𝑥0𝑍^02𝑑𝑐_𝑖∶𝑁=𝑋_{𝑖+1∶𝑁}−𝑋_𝑖∶𝑁1/𝑁𝑓(𝑋_𝑖∶𝑁). (3.12)In Theorem 1.5, we use the same scale factor for all differences in a bin (see (3.4)). As we assume the first and second derivatives of 𝑓 are uniformly bounded, (2.5) and (2.6) imply that for 𝑋_𝑖∶𝑁∈[𝑎_𝑘;𝑁,𝑏_𝑘;𝑁],𝑓(𝑋_𝑖∶𝑁)=𝑓(𝑎_𝑘;𝑁)+𝑂(𝑏_𝑘;𝑁−𝑎_𝑘;𝑁)=𝑓(𝑎_𝑘;𝑁𝑁)+𝑂(^𝛿−1𝑓(𝑎_𝑘;𝑁)+𝑁^2𝛿−2), (3.13)and the Big-Oh constants are independent of 𝑘. As we assume that 𝑓 satisfies (1.12), the error term is negligible.
Thus, our assumptions on 𝑓 imply that 𝑓 is basically constant on each bin, and we may replace the local rescaling factor 𝑓(𝑋_𝑖∶𝑁) with the bin rescaling factor 𝑓(𝑎_𝑘;𝑁). Thus, each bin of normalized differences has the same shift in its shifted exponential behavior. Therefore all the shifts reinforce, and the digits of all the normalized differences exhibit shifted exponential behavior as 𝑁→∞.

As an example of Theorem 1.7, in Figure 1 we consider 500,000 independent random variables drawn from the Pareto distribution with exponent𝑎=4+³19−3√33+³19+3√333 (3.14)(we chose 𝑎 to make the variance equal 1). We study the distribution of the digits of the differences in base 10. The amplitude is about .018, which is the amplitude of the shifted exponential behavior of Theorem 1.1 (see the equation in [23, Theorem 2] or (1.5) of Theorem 1.1).

Remark 3.5. The universal behavior of Theorem 1.7 suggests that if we are interested in the behavior of the digits of all the differences, the natural quantity to study is the normalized differences. For any distribution with uniformly bounded first and second derivatives and a second-order Taylor series expansion at each point, we obtain shifted exponential behavior.

Appendices

A. Proof of Theorem 1.1

To prove Theorem 1.1, it suffices to study the distribution of log_𝐵𝜁mod1 when 𝜁 has the standard exponential distribution (see (1.4)). We have the following useful chain of equalities. Let [𝑎,𝑏]⊂[0,1]. Then,Prob(log_𝐵𝜁mod1∈[𝑎,𝑏])=^∞_𝑘=−∞Prob(log_𝐵=𝜁∈[𝑎+𝑘,𝑏+𝑘])_{∞𝑘=−∞}Prob(𝜁∈[𝐵^𝑎+𝑘,𝐵^𝑏+𝑘=])_{∞𝑘=−∞}(𝑒^{−𝐵𝑎+𝑘}−𝑒^{−𝐵𝑏+𝑘}). (A.1)It suffices to investigate (A.1) in the special case when 𝑎=0, as the probability of any interval [𝛼,𝛽] can always be found by subtracting the probability of [0,𝛼] from [0,𝛽]. We are, therefore, led to study, for 𝑏∈[0,1], the cumulative distribution function of log_𝐵𝜁mod1,𝐹_𝐵(𝑏)∶=Prob(log_𝐵𝜁mod1∈[0,𝑏])=^∞_𝑘=−∞(𝑒^−𝐵𝑘−𝑒^{−𝐵𝑏+𝑘}). (A.2)This series expansion converges rapidly, and Benford behavior for 𝜁 is equivalent to the rapidly converging series in (A.2) equalling 𝑏 for all 𝑏.

As Benford behavior is equivalent to 𝐹_𝐵(𝑏) equals 𝑏 for all 𝑏∈[0,1], it is natural to compare 𝐹^′_𝐵(𝑏) to 1. If the derivatives were identically 1, then 𝐹_𝐵(𝑏) would equal 𝑏 plus some constant. However, (A.2) is zero when 𝑏=0, which implies that this constant would be zero. It is hard to analyze the infinite sum for 𝐹_𝐵(𝑏) directly. By studying the derivative 𝐹^′_𝐵(𝑏), we find a function with an easier Fourier transform than the Fourier transform of 𝑒^−𝐵𝑢−𝑒^{−𝐵𝑏+𝑢}, which we then analyze by applying Poisson summation.

We use the fact that the derivative of the infinite sum 𝐹_𝐵(𝑏) is the sum of the derivatives of the individual summands. This is justified by the rapid decay of the summands (see, e.g., [28, Corollary 7.3]). We find𝐹^′_𝐵(𝑏)=^∞_𝑘=−∞𝑒^{−𝐵𝑏+𝑘}𝐵^𝑏+𝑘log𝐵=^∞_𝑘=−∞𝑒^{−𝛽𝐵𝑘}𝛽𝐵^𝑘log𝐵, (A.3)where for 𝑏∈[0,1], we set 𝛽=𝐵^𝑏.

Let 𝐻(𝑡)=𝑒^{−𝛽𝐵𝑡}𝛽𝐵^𝑡log𝐵; note 𝛽≥1. As 𝐻(𝑡) is of rapid decay in 𝑡, we may apply Poisson summation (e.g., [29]). Thus,^∞_𝑘=−∞𝐻(𝑘)=^∞_𝑘=−∞𝐻(𝑘), (A.4)where 𝐻 is the Fourier transform of 𝐻: 𝐻(𝑢)=∫^∞_−∞𝐻(𝑡)𝑒^{−2𝜋𝑖𝑡𝑢}𝑑𝑡. Therefore,𝐹^′_𝐵(𝑏)=^∞_𝑘=−∞𝐻(𝑘)=^∞_𝑘=−∞𝐻(𝑘)=^∞_𝑘=−∞^∞_−∞𝑒^{−𝛽𝐵𝑡}𝛽𝐵^𝑡log𝐵⋅𝑒^{−2𝜋𝑖𝑡𝑘}𝑑𝑡. (A.5)Let us change variables by taking 𝑤=𝐵^𝑡. Thus, 𝑑𝑤=𝐵^𝑡log𝐵𝑑𝑡 or 𝑑𝑤/𝑤=log𝐵𝑑𝑡. As 𝑒^{−2𝜋𝑖𝑡𝑘}=(𝐵^𝑡/log𝐵)^{−2𝜋𝑖𝑘}=𝑤^{−2𝜋𝑖𝑘/log𝐵}, we have𝐹^′_𝐵(𝑏)=^∞_𝑘=−∞^∞₀𝑒^−𝛽𝑤𝛽𝑤⋅𝑤^{−2𝜋𝑖𝑘/log𝐵}𝑑𝑤𝑤=^∞_𝑘=−∞𝛽^{2𝜋𝑖𝑘/log𝐵∞}₀𝑒^−𝑢𝑢^{−2𝜋𝑖𝑘/log𝐵}=𝑑𝑢_{∞𝑘=−∞}𝛽^{2𝜋𝑖𝑘/log𝐵}Γ(1−2𝜋𝑖𝑘log𝐵), (A.6)where we have used the definition of the Γ-functionΓ(𝑠)=^∞₀𝑒^−𝑢𝑢^𝑠−1𝑑𝑢,Re(𝑠)>0. (A.7)As Γ(1)=1, we have𝐹^′_𝐵(𝑏)=1+^∞_𝑚=1[𝛽^{2𝜋𝑖𝑚/log𝐵}Γ(1−2𝜋𝑖𝑚log𝐵)+𝛽^{−2𝜋𝑖𝑚/log𝐵}Γ(1+2𝜋𝑖𝑚log𝐵)]. (A.8)

Remark A.1. The above series expansion is rapidly convergent, and shows the deviations of log_𝐵𝜁mod1 from being equidistributed as an infinite sum of special values of a standard function. As 𝛽=𝐵^𝑏, we have 𝛽^{2𝜋𝑖𝑚/log𝐵}=cos(2𝜋𝑚𝑏)+𝑖sin(2𝜋𝑚𝑏), which gives a Fourier series expansion for 𝐹^′(𝑏) with coefficients arising from special values of the Γ-function.

We can improve (A.8) by using additional properties of the Γ-function. If 𝑦∈ℝ, then from (A.7), we have Γ(1−𝑖𝑦)=Γ(1+𝑖𝑦) (where the bar denotes complex conjugation). Thus, the 𝑚th summand in (A.8) is the sum of a number and its complex conjugate, which is simply twice the real part. We have formulas for the absolute value of the Γ-function for large argument. We use (see [30, page 946, equation (8.332)]) that|Γ(1+𝑖𝑥)|²=𝜋𝑥=sinh(𝜋𝑥)2𝜋𝑥𝑒^𝜋𝑥−𝑒^−𝜋𝑥. (A.9)Writing the summands in (A.8) as 2Re(𝑒^{−2𝜋𝑖𝑚𝑏}Γ(1+2𝜋𝑖𝑚/log𝐵)), (A.8) becomes𝐹^′_𝐵(𝑏)=1+2^𝑀−1_𝑚=1Re(𝑒^{−2𝜋𝑖𝑚𝑏}Γ(1+2𝜋𝑖𝑚log𝐵))+2^∞_𝑚=𝑀Re(𝑒^{−2𝜋𝑖𝑚𝑏}Γ(1+2𝜋𝑖𝑚log𝐵)). (A.10)The rest of the claims of Theorem 1.1 follow from simple estimation, algebra, and trigonometry.

With constants as in the theorem, if we take 𝑀=1 and 𝐵=𝑒 (resp., 𝐵=10) the error is at most .00499 (resp., .378), while if 𝑀=2 and 𝐵=𝑒 (resp., 𝐵=10), the error is at most 3.16⋅10⁻⁷ (resp., .006). Thus, just one term is enough to get approximately five digits of accuracy base 𝑒, and two terms give three digits of accuracy base 10. For many bases, we have reduced the problem to evaluate Re(𝑒^{−2𝜋𝑖𝑏}Γ(1+2𝜋𝑖/log𝐵)). This example illustrates the power of Poisson summation, taking a slowly convergent series expansion and replacing it with a rapidly converging one. Corollary A.2. Let 𝜁 have the standard exponential distribution. There is no base 𝐵>1 such that 𝜁 is Benford base 𝐵.Proof. Consider the infinite series expansion in (1.5). As 𝑒^{−2𝜋𝑖𝑚𝑏} is a sum of a cosine and a sine term, (1.5) gives a rapidly convergent Fourier series expansion. If 𝜁 were Benford base 𝐵, then 𝐹^′_𝐵(𝑏) must be identically 1; however, Γ(1+2𝜋𝑖𝑚/log𝐵) is never zero for 𝑚 a positive integer because its modulus is nonzero (see (A.9)). As there is a unique rapidly convergent Fourier series equal to 1 (namely, 𝑔(𝑏)=1; see [29] for a proof), our 𝐹^′_𝐵(𝑏) cannot identically equal 1.

B. Analyzing 𝑁^1−𝛿 Intervals Simultaneously

We show why in addition to 𝜖>0 we also needed 𝜖>1/3−𝛿/2 when we analyzed 𝑁^1−𝛿 intervals simultaneously in (3.2); we thank one of the referees for providing this detailed argument.

Let 𝑌₁,…,𝑌_𝑁 be i.i.d.r.v. with 𝔼[𝑌_𝑖]=0, Var(𝑌_𝑖)=𝜎², 𝔼[|𝑌_𝑖|³]<∞, and set 𝑆_𝑁=(𝑌₁+⋯+𝑌_𝑁)/√𝑁𝜎². Let Φ(𝑥) denote the cumulative distribution function of the standard normal. Using a (nonuniform) sharpening of the Berry-Esséen estimate (e.g., [31]), we find that for some constant 𝑐>0,|Prob(𝑆_𝑁≤𝑥)−Φ(𝑥)|≤𝑐𝔼[|𝑌₁|³]𝜎³√𝑁(1+|𝑥|)³,𝑥∈ℝ,𝑁≥1. (B.1)Taking 𝑌_𝑖=𝑤_𝑖−𝑁^𝛿−1, where 𝑤_𝑖 is defined by (2.1), yields𝑆_𝑁=𝑀_𝑁−𝑁^𝛿𝑁^𝛿(1−𝑁^𝛿−1),𝜎²𝑁^𝛿−1(1−𝑁^𝛿−1),𝔼[|𝑌_𝑖|³]≤2𝑁^𝛿−1. (B.2)Thus, (B.1) becomes𝑀|Prob(^𝑁−𝑁^𝛿𝑁^𝛿(1−𝑁^𝛿−1)≤𝑥)−Φ(𝑥)|≤3𝑐𝑁^−𝛿/2(1+|𝑥|)³, (B.3)for all 𝑁≥𝑁₀ (for some 𝑁₀ sufficiently large, depending on 𝛿).

For each 𝑁, 𝑘, and 𝜖 consider the event𝐴_{𝑁,𝑘,𝜖}𝑀={^𝑘;𝑁−𝑁^𝛿𝑁^𝛿(1−𝑁^𝛿−1)∈[−𝑁^𝜖,𝑁^𝜖]}. (B.4)Then, as 𝑁→∞, we haveProb(^𝑁1−𝛿_𝑘=1𝐴_{𝑁,𝑘,𝜖})→1, (B.5)provided that^𝑁1−𝛿_𝑘=1Prob(𝐴^𝑐_{𝑁,𝑘,𝜖})→0, (B.6)as 𝑁→∞. Using (B.3) givesProb(𝐴^𝑐_{𝑁,𝑘,𝜖})≤6𝑐𝑁^−𝛿/2(1+𝑁^𝜖)³+2(1−Φ(𝑁^𝜖))≤6𝑐𝑁^{−𝛿/2−3𝜖}+2𝜋𝑁^−𝜖𝑁exp(−^2𝜖2) (B.7)(e.g., [32]). Thus, the sum in (B.6) is at most6𝑐𝑁^{1−3𝛿/2−3𝜖}+2𝜋𝑁^{1−𝛿−𝜖}𝑁exp(−^2𝜖2), (B.8)and this is 𝑂(1) provided that 𝜖>0 and 𝜖>1/3−𝛿/2.

Acknowledgments

The authors would like to thank Ted Hill, Christoph Leuenberger, Daniel Stone, and the referees for numerous helpful comments. S. J. Miller was partially supported by NSF (Grant no. DMS-0600848).