The Newcomb-Benford law: Scale invariance and a simple Markov process based on it

The Newcomb-Benford law, also known as the first-digit law, gives the probability distribution associated with the first digit of a dataset, so that, for example, the first significant digit has a probability of $30.1$ % of being $1$ and $4.58$ % of being $9$. This law can be extended to the second and next significant digits. This article presents an introduction to the discovery of the law, its derivation from the scale invariance property, as well as some applications and examples, are presented. Additionally, a simple model of a Markov process inspired by scale invariance is proposed. Within this model, it is proved that the probability distribution irreversibly converges to the Newcomb-Benford law, in analogy to the irreversible evolution toward equilibrium of physical systems in thermodynamics and statistical mechanics.


I. INTRODUCTION
In the late 19th century, an astronomer and mathematician visits his institution's library and consults a table of logarithms to perform certain astronomical calculations. As on previous occasions, he is struck by the fact that the first pages (those corresponding to numbers that start at 1) are much more worn than the last ones (corresponding to numbers that start at 9). Intrigued, this time he decides not to overlook the matter. He closes his eyes to concentrate, sketches a few calculations on a piece of paper, and finally smiles. He has found the answer and it turns out to be enormously simple and elegant.
A little over half a century later, a physicist and electrical engineer who is unaware of his predecessor's discovery, observes the same curious phenomenon on the pages of logarithm tables, and arrives at the same conclusion. Both have understood that, in a long list of records {r n } obtained from measurements or observations, the fraction p d of records beginning with the significant digit d = 1, 2, . . . , 9 is not p d = 1/9, as one might naively expect, but rather follows a logarithmic law. More specifically, p d = log 10 1 + 1 d , d = 1, 2, . . . , 9. (1) The numeric values of p d are shown in the second column of Table I. We see that the records that start with 1, 2, or 3 account for around 60 % of the total, while the other six digits must settle for the remaining 40 %. Our 19th century character is Simon Newcomb (Fig.  1) and he published his discovery in a modest two-page note. 1 The second character is Frank Benford (Fig. 2) and he wrote a 22-page article 2 in which, in addition to mathematically justifying Eq. (1), he showed its validity in the analysis of more than 20 000 first digits a) Electronic mail: anburgosr@alumnos.unex.es b) Electronic mail: andres@unex.es taken from sources as varied as river areas, populations of American cities, physical constants, atomic and molecular weights, specific heats, numbers taken from newspapers or Reader's Digest, postal addresses, . . . , and the series n −1 , √ n, n 2 , or n!, among others, with n ∈ 1, 100 . With such an overwhelming evidence, it is not surprising that Eq. (1) is usually known as Benford's law (or first-digit law), even though it was found nearly sixty years earlier by Newcomb. This is one more manifestation of the so-called Stigler's law, according to which no scientific discovery is named after the person who discovered it in the first place. In fact, as Stigler himself points out, 3 the law that bears his name was actually spelled out in a similar way twenty-three years earlier by the American sociologist Robert K. Merton. In order not to fall completely into Stigler's law, many authors refer to Eq. (1) as Newcomb-Benford's law and that is the convention (by means of the acronym NBL) that we will follow in this article.
While the literature on the NBL in specialized journals is vast, 4 and several books on the topic are available, 5,6 the number of papers in general or science education journals is much scarcer. [7][8][9] In this paper, apart from revis- FIG. 1. Simon Newcomb (1835-1909. FIG. 2. Frank Benford (1883-1948 iting the NBL and illustrating it with a few examples, we construct a Markov-chain model inspired by the invariance property of the NBL under the operation of a change of scale. The remainder of the paper is organized as follows. The connection between the NBL (and some of its generalizations) with the property of scale invariance is formulated in Sec. II. This is followed by a few examples in Sec. III. The most original part is contained in Sec. IV, where our Markov-chain model is proposed and solved. Finally, the paper is closed in Sec. V with some concluding remarks. For the interested reader, Appendices A-C contain the most technical and mathematical parts of Secs. II and IV.

II. ORIGIN OF THE LAW
Often times, when one first speaks to a friend, relative, or even a colleague about the NBL, their first reaction is usually of skepticism. Why is the first digit not evenly distributed among the nine possible values? A simple argument shows that, if a robust distribution law exists, it cannot be the uniform distribution whatsoever.
Imagine a long list of river lengths, mountain heights, and country surfaces, for example. It is possible that the lengths of the rivers are in km, the heights of the mountains in m, and the surfaces of countries in km 2 , but they could also be expressed in miles, feet, or acres, respectively. Will the distribution p d depend on whether we use some units or others, or even if we mix them? It seems logical that the answer should be negative; that is, the p d distribution should be (statistically) independent of the units chosen; in other words, it is expected that the p d distribution is invariant under a change of scale. The uniform distribution p d = 1 9 obviously does not follow that property of invariance. Suppose we start from a dataset in which all the values of the first digit are equally represented. If we multiply all the records in the dataset by 2, we can see that those records that started before with 1, 2, 3, and 4 then start with either 2 or 3, either 4 or 5, either 6 or 7, and either 8 or 9, respectively. On the other hand, all those that started with 5, 6, 7, 8, or 9 start now with 1. All those possibilities are schematically shown in Fig. 3. Therefore, if p d = 1 9 initially, then p 1 = 5 9 and p 2 + p 3 = p 4 + p 5 = p 6 + p 7 = p 8 + p 9 = 1 9 after doubling all the records, thus destroying the initial uniformity. Thus, the most identifying hallmark of the NBL is that it must be applied to records that have units or, as Newcomb himself writes, 1 "As natural numbers occur in nature, they are to be considered as the ratios of quantities." While relatively formal mathematical proofs of the NBL can be found in the literature, [10][11][12] in Appendix A, we present a sketch of a simpler, physicist-style derivation of the law by imposing invariance under a change of scale. 13 Equation (1) can be generalized beyond the first digit. The probability that the m first digits of a record match the ordered string (d 1 , d 2 , . . . , d m ), where d 1 ∈ {1, 2, . . . , 9} and d i ∈ {0, 1, 2, . . . , 9} if i ≥ 2, is given by (see Appendix A) p d1,d2,...,dm = log 10 As an example, the probability that the first three digits of a record form precisely the string (3,1,4) is p 3,1,4 = log 10 (1 + 1/314) = 0.001 38.
Once we have p d1,d2,...,dm , we can calculate the probability p In Table I, the law for the first digit, p d , is accompanied by the laws for the second, third, and fourth digits, obtained from Eqs. (2) and (3). As the digit moves to the right, the probability distribution becomes less and less disparate. In Fig. 3 we saw that, when multiplying a dataset {r n } by 2, part of the records that started with d = 1, 2, 3, 4, specifically those that start with the first two digits (d, 0), (d, 1), (d, 2), (d, 3), or (d, 4), will start with 2d, while the rest will start with 2d + 1. Let us call α d the fraction of records that, initially starting with d = 1, 2, 3, 4, start with 2d by doubling all the records. Thus, If the dataset fulfills the NBL, then one has α 1 = log 10 3 2 log 10 2 ≃ 0.584 96, α 2 = log 10 5 4 We will use these values later in Sec. IV.

III. APPLICATIONS AND EXAMPLES
The applications and verifications of the NBL are numerous and cover topics as varied and prosaic as the study of the genome, 14 atomic, 15 nuclear, 8,16 and particle 17 physics, astrophysics, 18 quantum correlations, 19 toxic emissions, 20 biophysics, 21 medicine, 22 dynamical systems, 23 distinction of chaos from noise, 24 statistical physics, 25 scientific citations, 26 tax audits, 5,27 electoral 28 or scientific 29 frauds, gross domestic product, 30 stock market, 9,31 inflation data, 32 climate change, 33 world wide web, 34 internet traffic, 35 social networks, 36 textbook exercises, 37 image processing, 38 religious activities, 39 dates of birth, 40 hydrology and geology, 41 fragmentation processes, 42 the first letters of words, 43 or even  Other examples can be seen in the link of Ref. 45. In this section, we will present some additional examples.
Let us start with one of the situations that Benford himself studied in his classic paper: 2 city populations. Using data from the Spanish National Institute of Statistics, 46 we have considered the 2019 population of the 165 municipalities in the province of Badajoz (plus the total population of the province of Badajoz), of the 223 municipalities of the province of Cáceres (plus the total population of the province of Cáceres), and the total population of the 388 municipalities of the Spanish region of Extremadura, which encompasses the provinces of Badajoz and Cáceres, (plus the total populations of the provinces of Badajoz and Cáceres). We have also considered the population (according to the 2016 census) of the 8 110 Spanish municipalities. For all these lists of records, we have analyzed the frequency of those starting with d = 1, 2, . . . 9 and the results are compared in Fig. 4. There is a good general agreement between the population data and the NBL predictions. This is interesting, since it is not obvious that the distribution of the number of inhabitants of municipalities should be invariant under a change of scale.
Let us now turn to two examples from astronomy. In the first one, we take the distance from Earth (in lightyears and in parsecs) to the 300 brightest stars. 47 In the second case, the data are the daily number of sunspots from 1818 to the present. 48 As seen in Fig. 5, the distances between our planet and the 300 brightest stars agree reasonably well with the NBL, despite the fact that the list is not excessively long (only 300 data points) and that there are "local" deviations (for example, p 6 < p 7 in the two choices of units). This general agreement was to be expected, since the distribution of digits in distances (which are expressed in units) is a clear example of invariance under a change of scale. However, in the case of the daily number of sunspots, quantitative (although not qualitative) differences are observed with the NBL, especially in the d = 1, 3, 4, and 5 cases. It should be noted that, although the series is rather long (more than 59 000 records, after excluding days without data or with zero spots), each record only has one, two, or three digits (the maximum number of sunspots was 528 and corresponded to August 26, 1870), thus producing a certain bias to records starting with 1. Moreover, the daily number of sunspots cannot be expressed in units and therefore may not be scale-invariant. Finally, we have analyzed the prices of 1 016 items from a chain of fashion retailers 49 and of 1 373 products from a chain of hypermarkets 50 . The results are shown in Fig. 6. In this case, the discrepancies with the NBL are more pronounced. Although the highest frequencies occur for d = 1 and d = 2, the observed values of p d do not decrease monotonically with increasing d.
In the case of the fashion retailers, we have p 4 > p 3 and p 9 > p 8 > p 6 > p 7 ; in the prices of the chain of hypermarkets, p 8 > p 9 > p 7 > p 6 . In principle, one might think that, since they can be expressed in euro, pound, dollar, peso, ruble, yen, . . . , prices should satisfy the property of invariance under a change of scale inherent to the NBL. However, commercial and artificial pricing strategies must be superimposed on this invari- ance, which generates relevant deviations with respect to the NBL.

IV. A SIMPLE MODEL OF A MARKOV CHAIN BASED ON THE SCALE INVARIANCE PROPERTY OF THE NEWCOMB-BENFORD DISTRIBUTION
As already said, NBL, Eq. (1), is invariant under a change of scale; that is, if we start from a dataset {r n } that fulfills the NBL and multiply all the records by a constant λ, the resulting dataset {λr n } still fulfills the NBL.
At this point, we construct a simplified dynamical model in which the four unknown and time-dependent parameters α d in Eqs. (6) are replaced by fixed constants. In matrix notation,  where p(t) = (p 1 (t), p 2 (t), . . . , p 9 (t)) † is a column vector  (7). The dashed line is proportional to |a2,3| 2t (see Appendix C).

( † denotes transposition) and
is a fixed square matrix. Equation (7) fits the canonical description of Markov chains, 52 where W is the socalled transition matrix, and {α d } correspond to transition probabilities. In this way, we may forget about the meaning of {p d (t)} as the first-digit distribution of the dataset {r n (t)} and look at the numerals 1, 2, . . . , 9 as labels of nine distinct internal states of a certain physical system which follow each other in the prescribed order sketched by Fig. 3.
Note that the matrix W is singular, that is, not invertible. This implies the irreversible character of the transition {p d (t)} → {p d (t + 1)}. The stationary solution p * to Eq. (7) satisfies the condition p * = W · p * . Such a stationary solution will be an attractor of our Markov process if lim t→∞ p(t) = p * for any initial condition p(0). This property, along with the exact solution of Eq. (7), is proved in Appendix B. If we choose for α d the values given by Eqs. (5), the stationary solution coincides exactly with the NBL, as could be expected. This will be the choice made in the rest of this section.
To illustrate the irreversible evolution of p(t) toward p * , we are going to consider two different initial conditions. First, we start from a uniform distribution, that is, p d (0) = 1 9 . The result is shown in Fig. 8, where we see that the evolution is oscillatory (see Appendix B for an explanation) and the oscillations are rapidly damped to attain the stationary solution. As a second example, we take an NBL inverted initial distribution, that is, p d (0) = p * 10−d , so that, initially, state 9 is the most probable and state 1 is the least probable. In this second example, as shown in Fig. 9, the initial oscillations are of greater amplitude but, as before, the stationary distribution is practically reached after a few iterations. Comparison between Figs. 7 and 8 shows that the main difference between our Markov model and the non-Markovian evolution of {p d (t)} in a real dataset is that the oscillation amplitudes are attenuated in the model and not in the original framework.
It seems convenient to characterize the evolution of the set of probabilities {p d (t)} obeying the Markov process [Eq. (7)] toward the attractor {p * d } by means of a single parameter that, in addition, evolves monotonically, thus representing the irreversibility of evolution. It is expected that these properties are verified by the Kullback-Leibler divergence, 53 which in our case is defined as This quantity represents the opposite of the Shannon entropy 54 of {p d (t)}, except that it is measured with respect to the stationary distribution {p * d }. In the context of our Markov-chain model, D KL is related to information theory. On the other hand, the mathematical structure of D KL can be extended to physical systems, thus providing a statistical-mechanical basis to the thermodynamic entropy, 54,55 as exemplified below. Figure 10 shows the evolution of D KL (t) for the same initial conditions as in Figs. 8 and 9. Both cases confirm the desired monotonic evolution of D KL (t). Also, the asymptotic decay D KL (t) → 0 occurs essentially exponentially with a rate independent of the initial probability distribution. This is proved in Appendix C, as well as the monotonicity property with the equality not holding for two successive times unless p d (t) = p * d , in which case D KL = 0. Thus, we can say that the NBL in our Markov model plays a role analogous to the equilibrium state in thermodynamics and statistical mechanics. In this analogy, the "entropy" of the out-of-equilibrium system would be (except for a constant) S = −D KL , so that S increases irreversibly in the evolution toward equilibrium.
This analogy is put forward in Table II, where we compare a system described by the Markov chain [Eq. (7)] to a classical dilute gas as an example of a well-known physical system. In both cases, a statistical description is constructed by introducing the adequate probability TABLE II. Analogy between a classical dilute gas (as a prototypical physical system) and a system described by the Markov chain, Eq. (7). In the expression of the Maxwell-Boltzmann distribution fMB( v), m is the mass of a molecule, T is the gas temperature, and kB is the Boltzmann constant. In the Boltzmann equation, J[ v|f (t), f (t)] is the collision operator.

Dilute gas Markov chain Probability distribution
Velocity distribution function: f ( v, t) Probability of the internal state d: p d (t) Normalization distribution: the velocity distribution function (gas) and the probability of the internal state d (Markov chain); while f ( v, t) is continuous in both v and t, p d (t) is discrete in d and t. The evolution equation for the probability distribution function is the Boltzmann equation (gas, under the molecular chaos ansatz 56 ) and Eq. (7) (Markov chain). The equilibrium state is characterized by the Maxwell-Boltzmann distribution 57 f MB ( v) (gas) and the NBL distribution p * d (Markov chain). In both classes of systems the nonequilibrium entropy functional S(t) can be defined, within a constant, as the opposite of the Kullback-Leibler divergence 53,58 of the equilibrium distribution from the time-dependent one. Boltzmann's celebrated H-theorem 56,59 (gas) and the result presented in Eq. (10) (Markov chain) show that the entropy S(t) never decreases, irreversibly evolving toward its equilibrium value.

V. CONCLUDING REMARKS
One of the main goals of this article was to show that, contrary to what might be initially thought, the first significant digit d of a dataset extracted from measurements or observations of the real world is not evenly distributed among the nine possible values, but typically the frequency is higher for d = 1 and decreases as d increases. The NBL, Eq. (1), gives a mathematical expression to this empirical fact, although it does not always need to be rigorously verified due to statistical fluctuations (as happens with populations in Fig. 4 and with distances in Fig. 5), limitations of the sample (as in the sunspot case in Fig. 5), or an artificial bias (as in the case of prices of articles in Fig. 6). It is to be expected that, except for unavoidable statistical fluctuations, the law is fulfilled in datasets in which quantities are measured in units, so that the distribution of the first digit is independent of the units chosen due to its invariance under a change of scale. More generally, the NBL is satisfied when the mantissa of the logarithms of the considered data is uniformly distributed. That makes lists not directly related to physical quantities, such as Fibonacci numbers or powers of 2, also satisfy the NBL.
Gaining inspiration from the scale invariance property of the NBL, we have constructed a Markov-chain model for a nine-state system that evolves in time according to the scheme depicted in Fig. 3. The initialvalue model can be exactly solved, the solution showing an irreversible relaxation toward the NBL probability distribution. Moreover, we have proved that the associated (relative) entropy monotonically increases, in analogy with the second law of thermodynamics in physical systems.
Until the 1970s (which is when scientific pocket calculators began to be used), physicists used tables of logarithms (or their application in slide rules) for small everyday scientific calculations. Those calculations are nowadays performed on pocket calculators, cellular phones, or personal computers with a wide variety of existing mathematical programs. Since the data that are manipulated in physics are extracted from "real" situations, such as experiments, models, physical constants, . . . , we can conclude, as a tribute to Newcomb and Benford and their logarithmic tables, that the keyboard button 1 will be the one that presents the greatest wear and tear, while that of 9 will be the least used. Consider a long list of records {r n } that, without loss of generality for the matter at hand, we will assume positive. Each record can be written in the form r n = x n × 10 kn , where k n is an integer and x n ∈ [1, 10) is the significand. Obviously, it is the distribution of the significand that is relevant for the NBL. The significand x n is directly related to the mantissa µ n of the decimal logarithm of r n , that is, log 10 r n = k n + µ n , where µ n = log 10 x n ∈ [0, 1).
Let P x (x)dx be the probability that the significand lies between x and x + dx, so that the normalization condition is 10 1 dx P x (x) = 1. The probability that the first significant digit of any record is d is then given by the integral More generally, if we consider an ordered string (d 1 , d 2 , . . . , d m ) made up of the first m digits, where d 1 ∈ {1, 2, . . . , 9} and d i ∈ {0, 1, 2, . . . , 9} if i ≥ 2, it is obvious that the records whose m first digits match the string (d 1 , d 2 , . . . , d m ) will be those whose significand x is greater than or equal to d 1 + d 2 × 10 −1 + · · · + d m × 10 −(m−1) and less than d 1 + d 2 × 10 −1 + · · · + (d m + 1) × 10 −(m−1) . Consequently, If the distribution P x (x) is actually invariant under a change of scale, that means that P x (λx) = f (λ)P x (x) with arbitrary λ. Taking into account the normalization condition in the form 10λ λ d(λx) P x (λx) = 1, it follows that necessarily f (λ) = λ −1 , that is, P x (λx) = λ −1 P x (x). Differentiating both sides of the equation with respect to λ and then taking λ = 1, we easily obtain xP ′ x (x) = −P x (x), which, according to Euler's theorem, implies that P x (x) is a homogeneous function of degree −1, that is, P x (x) ∝ x −1 . The constant of proportionality is obtained from the normalization condition, finally yielding This is the unique distribution of significands that is invariant under a change of scale. From Eq. (A3), and applying Eqs. (A1) and (A2), it is straightforward to obtain NBL, Eq. (1), and its generalization, Eq. (2). As a consistency test of Eq. (2), it is easy to check that p d1,d2,...,dm−1 = 9 dm=0 p d1,d2,...,dm = log 10 It is interesting to note that the inverse law for the significand, Eq. (A3), implies a uniform law for the mantissa (and vice versa). Let P µ (µ)dµ be the probability that the mantissa lies between µ and µ + dµ. Since P µ (µ)dµ = P x (x)dx and dµ = (x −1 / ln 10)dx, Eq. (A3) gives us P µ (µ) = 1. In Newcomb's words, 1 "The law of probability of the occurrence of numbers is such that all mantissae of their logarithms are equally probable." An immediate consequence is that, if µ is a random variable uniformly distributed between 0 and 1, then the random variable x = 10 µ fulfills the NBL. This is in fact an easy way to generate a list of random records matching the NBL.
There are deterministic series that also satisfy the NBL. Suppose the series {r n = aα n + bβ n , n = 1, 2, . . .}, where a, b, α, and β are real numbers with a > 0, α > β ≥ 0, and log 10 α = irrational. 12 In that case, lim n→∞ log 10 r n = n log 10 α + log 10 a has a uniformly distributed mantissa, so {r n } satisfies the NBL. That includes, for example, the series {2 n }, {3 n }, and {F n }, where F n = 1 √ 5 ϕ n − (−ϕ −1 ) n are the Fibonacci numbers, ϕ ≡ 1 2 1 + √ 5 being the golden ratio. Similarly, the series {n!} also verifies the law. 12 Another important property is that if a list {r n } fulfills the NBL, so does the list {r a n }, a being a real number. Indeed, if log 10 r n = k n + µ n , the mantissa µ n being uniformly distributed, then the mantissa of log 10 r a n = a(k n + µ n ) is also evenly distributed. This is directly related to the fact that the NBL is not only invariant under a change of scale but also under base change, 11 as would be expected, given the artificial character of the decimal base. To see it, let us assume a base b and build the list {r a n }, with a = log 10 b, from a list {r n } that fulfills the NBL. In that case, r a n = y n × b kn , where y n = x a n ∈ [1, b) is the significand of r a n in the base b. The probability distribution P y (y) is related to the distribution P x (x) through the equation P y (y)dy = P x (x)dx, so that Eq. (A3) leads to Therefore, NBL in an arbitrary base b takes the form Therefore, only eight of the probabilities {p d , d = 1, 2, . . . , 9} are linearly independent, so we can eliminate one of them. If we choose p 9 = 1 − 8 d=1 p d , Eq. (6a) gives us p 1 (t + 1) = 1 − p 1 (t) − p 2 (t) − p 3 (t) − p 4 (t). Although it is not strictly necessary, it is mathematically convenient to split the column vector p(t) into the subsets p I (t) ≡ (p 1 (t), p 2 (t), p 3 (t), p 4 (t)) † , p II (t) ≡ (p 5 (t), p 6 (t), p 7 (t), p 8 (t)) † , and p 9 (t). As a consequence, the matrix equation (7) can be decomposed into In this way, one deals with two independent 4×4 matrices instead of the 9 × 9 transition matrix introduced in Eq. (8). Moreover, only the equation for the projected vector p I in Eq. (B1) needs to be solved. Once the solution is obtained, the solution for the complementary projected vector p II is straightforward. By recursion, it is easy to check that the solution to the initial-value problem posed by Eq. (B1) is p II (t) = B · (I − A t−1 ) · p * I + B · A t−1 · p I (0), (B3b) where I is the identity matrix and p * II = B · p * I = is the unique stationary solution. From the normalization condition, one has p * 9 = α 1 α 2 (1 − α 4 )/(3 + α 1 α 2 ). Note that, as seen from Eqs. (B3), the initial values p II (0) do not influence the evolution of either p I (t) or p II (t).
The stationary solution will be an attractor if lim t→∞ p I (t) = p * I and lim t→∞ p II (t) = p * II for any initial condition p I (0). According to Eqs. (B3), this implies lim t→∞ A t = 0.