key: cord-0307967-5b6wetix authors: Burgos, Andrea; Santos, Andr'es title: The Newcomb-Benford law: Scale invariance and a simple Markov process based on it date: 2021-01-28 journal: nan DOI: 10.1119/10.0004957 sha: c7c649c924ea588d5975d29112ed9eee96bd21d4 doc_id: 307967 cord_uid: 5b6wetix 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. 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. 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. 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. 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. It is tempting to conjecture that the NBL should be an attractor of this process. This would mean that, if we started from an initial dataset {r n (0)} that does not fulfill the NBL and generated new sets {r n (t)} = {λ t r n (0)} at times t = 1, 2, . . . by multiplying each successive dataset by λ (other than a fractional power of 10), the first-digit distribution {p d (t)} of the generated sets would converge toward the NBL. However, this is not the case. As illustrated in Fig. 7 for the case λ = 2 and an initial dataset of 10 4 random numbers, the timedependent distribution {p d (t)} exhibits a quasiperiodic oscillatory pattern around the NBL distribution without any apparent amplitude attenuation. In fact, since 2 10 = 1 024 ≃ 10 3 , the distribution nearly returns to the uniform initial one at times t = 10, 20, 30, . . .. This behavior is reminiscent of the Poincaré recurrence time in mechanical systems and the associated Zermelo paradox about irreversibility. 51 In the transformation {r n (t)} → {r n (t + 1) = 2r n (t)}, the first-digit distribution changes, according to Fig. 3 , as p 1 (t + 1) = p 5 (t) + p 6 (t) + p 7 (t) + p 8 (t) + p 9 (t), (6a) p 2 (t + 1) = α 1 p 1 (t), p 3 (t + 1) = (1 − α 1 )p 1 (t), (6b) p 4 (t + 1) = α 2 p 2 (t), p 5 (t + 1) = (1 − α 2 )p 2 (t), (6c) where the fractions α d (d = 1, 2, 3, 4) are defined by Eq. (4) . Note that, in general, the fractions α d depend on the distributions of the first digit, p d (t), and of the first two digits, p d,d2 (t), of the dataset {r n (t)} [see Eq. (4)]. As a consequence, (i) Eqs. (6) do not make a closed set and (ii) the parameters α d depend on t. 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). 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. 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. 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. To check the above condition, let us obtain the four eigenvalues {a i , i = 0, 1, 2, 3} of A. It is easy to see that the characteristic equation is a(α 1 α 2 + a + a 2 + a 3 ) = 0. Therefore, the eigenvalues are a 0 = 0 and where ı is the imaginary unit and β ≡ 3 2 7 3 − 9α 1 α 2 + 9 − 42α 1 α 2 + 81α 2 1 α 2 2 1/3 . A t = U · D t · U −1 , t = 1, 2, . . . , (1−α1)a2 α1α2 (1−α1)a3 α1α2 From Eqs. (B5) it can be verified that |a 1 | < |a 2,3 | < 1 for 0 < α 1 α 2 < 1, so that lim t→∞ a t i = 0 for i = 1, 2, 3. Therefore, lim t→∞ D t = 0 and hence lim t→∞ A t = 0. This proves the character of the stationary distribution p * as an attractor of the dynamics. Moreover, since both the real eigenvalue (a 1 ) and the real part of the complex eigenvalues (a 2,3 ) are negative, the evolution of each p d (t) to p * d is oscillatory, as can be observed in Figs. 8 and 9. Note that, in the analysis above, the eigenvalues 0 (double), 1 − α 3 , and α 4 of the matrix B are not needed. If we choose α d = 1 2 , then Eqs. (B4) provide the stationary solution p * 1 = 4 13 ≃ 0.308, p 2 = p 3 = 2 13 ≃ 0.154, p 4 = p 5 = p 6 = p 7 = 1 13 ≃ 0.077, p 8 = p 9 = 1 26 ≃ 0.038. These values are not too different from those of the true NBL, as can be seen from Table I . On the other hand, the most natural choice is provided by Eqs. (5) , in which case the stationary solution coincides exactly with the NBL. The law of the anomalous numbers Stigler's law of eponymy Benford online bibliography Benford's Law: Applications for Forensic Accounting, Auditing, and Fraud Detection An Introduction to Benford's Law Benford's Law: Theory and Applications Benford's Law: Theory, The General Law Of Relative Quantities, And Forensic Fraud Detection Applications Forensic Analytics: Methods and Techniques for Forensic Accounting Investigations Significant figures of numbers in statistical tables Distribution of numbers and distribution of significant figures On the numbers of things and the distribution of first digits Benford's law and physical constants: The distribution of initial digits La misteriosa ley del primer dígito How do numbers begin? A simple explanation of Benford's law Benford's law: A "sleeping beauty" sleeping in the dirty pages of logarithmic tables Thermodynamics of Benford's first digit law An illustration of Benford's first digit law using alpha decay half lives The first digit phenomenon: A century-old observation about an unexpected pattern in many numerical tables applies to the stock market, census statistics and accounting data On the distribution of first significant digits On the probability that a random integer has initial digit A On the distribution of first significant figures The distribution of leading digits and uniform distribution mod 1 On Benford's law A statistical derivation for the significant-digit law Survival distributions satisfying Benford's law First digit law from Laplace transform A proof of first digit law from Laplace transform Uniform distribution, Benford's law and scale-invariance The mathematics of Benford's law: a primer Base-invariance implies Benford's law A basic theory of Benford's law Benford's law Making sense of microarray data distributions Benford's law and complex atomic spectra Benford's law and half-lives of unstable nuclei Benford's law and β-decay half-lives First digit distribution of hadron full width Do τ lepton branching fractions obey Benford's law? Newcomb-Benford law in astrophysical sources Empirical mantissa distributions of pulsars Benford's law in astronomy Hubble's law implies Benford's law for distances to galaxies Benford's distribution in extrasolar world: Do the exoplanets follow Benford's distribution? Benford's law in the Gaia universe Statistics of leading digits leads to unification of quantum correlations Benford analysis of quantum critical phenomena: First digit provides high finite-size scaling exponent while first two and further are not much better Quantum discord and its allies: a review of recent progress Assessing the accuracy of selfreported data: An evaluation of the toxics release inventory On the validation of the Newcomb-Benford law and the Weibull distribution in neuromuscular transmission Using the Benford's law as a first step to assess the quality of the cancer registry data Do dynamical systems follow Benford's law? Stochastic aspects of one-dimensional discrete dynamical systems: Benford's law One-dimensional dynamical systems and Benford's law Beyond Benford's law: Distinguishing noise from chaos The significant digit law in statistical physics Benford's law and citations, articles and impact factors of scientific journals Benford's law and articles of scientific journals: comparison of JCR and Scopus data Data diagnostics using secondorder tests of Benford's law Breaking the (Benford) law. Statistical fraud detection in campaign finance A firstdigit anomaly in the 2009 Iranian presidential election Statistical analysis of Brazilian electoral campaigns via Benford's law Not the first digit! Using Benford's law to detect fraudulent scientific data Benford's law first significant digit and distribution distances for testing the reliability of financial reports in developing countries Evidence for gross domestic product growth time delay dependence over foreign direct investment. A time-lag dependent correlation study Explaining the uneven distribution of numbers in nature: the laws of Benford and Zipf Tampering with inflation data: A Benford law-based analysis of national statistics in Argentina Technological improvements or climate change? Bayesian modeling of time-varying conformance to Benford's law Frequency of occurrence of numbers in the World Wide Web Benford's law behavior of internet traffic Benford's law applies to online social networks Benford's law: Textbook exercises and multiple-choice testbanks Images and Benford's law The Benford law behavior of the religious activity data Econophysics of a religious cult: The Antoinists in Belgium Breakdown of Benford's law for birth data Benford's law applied to hydrology data-Results and relevance to other geophysical data Benford's law in the natural sciences Applying Benford's law to volcanology Long-range properties and data validity for hydrogeological time series: The case of the Paglia river On the Newcomb-Benford law Benford's law and continuous dependent random variables Benford's law and first letter of words COVID-19, flattening the curve, and Benford's law On the authenticity of COVID-19 case figures The brightest stars Zermelo, Boltzmann, and the recurrence paradox Mechanical model for 2-state Markov chains On information and sufficiency Entropy, information, and computation The Second Law Reduced to Plain Common Sense The Mathematical Theory of Non-Uniform Gases The Boltzmann Equation and Its Applications Kinetic Theory of Gases in Shear Flows: Nonlinear Transport, Fundamental Theories of Physics Direct simulation for a homogeneous gas Thermodynamic time asymmetry and the Boltzmann equation Equilibrium solution of the Boltzmann equation A different proof of the Maxwell-Boltzmann distribution Derivation of the Maxwellian distribution from the microcanonical ensemble On deriving the Maxwellian velocity distribution Great moments in kinetic theory: 150 years of Maxwell's (other) equations Alternative derivation of the Maxwell distribution of speeds Time evolution of entropy, in various scenarios Boltzmann's H theorem applied to simulations of polymer interchange reactions Boltzmann's H-theorem and the assumption of molecular chaos Origin of the thermodynamic time arrow demonstrated in a realistic statistical system Appendix B: Exact solution of the Markov-chain model Note first that Eqs. (6) verify the normalization condition 9 d=1 p d (t + 1) = Appendix C: Properties of the Kullback-Leibler divergence in the Markov model The Kullback-Leibler divergence is defined by Eq. (9) . In order to analyze its asymptotic decay, let us consider times that are long enough so that the deviations δp d (t) ≡ p d (t) − p * d can be considered small. In this regime, we can expand Eq. (9) in a power series and retain the dominant terms. The result isOn the other hand, for times long enough, |a 1 | t ≪ |a 2,3 | t (note that |a 1 | = 0.4261 and |a 2,3 | = 0.8692), so that,This asymptotic behavior is represented in Fig. 10 .Finally, let us prove Eq. (10). According to Eq. (9),Evolution equations (6) and the stationarity condition p * = W · p * can be rewritten asInserting Eqs. (C4) into Eq. (C3) one gets(C5) Therefore,Given the stationary values {p * d , d = 5, . . . , 9}, the difference ∆D KL (t + 1) is a function of the 5 probabilities {p d (t), d = 5, . . . , 9}. To find the maximum value of ∆D KL (t + 1), we take the derivative. The unique solution to the extremum conditions 5, . . . , 9) , where 0 < γ < 1/ 9 d=5 p * d is arbitrary. In such a case, ∆D KL (t + 1) = 0. To see that this is actually a maximum value, suppose, for instance, that p d (t) = 0 except if d = d 0 (with d 0 = 5, . . . , 9) . In that case, it is easy to find ∆D KL (t + 1) = p d0 (t) ln p * d0 / 9 d=5 p * d < 0. We then conclude that ∆D KL (t + 1) ≤ 0, i.e., we obtain Eq. (10), the equality holding only if p d (t) = γp * d (d = 5, . . . , 9) . Note that, even though ∆D KL (t + 1) = 0 if p d (t)/p * d = constant for d = 5, . . . , 9, this proportionality property is broken down at t + 1 unless γ = 1. As a result, ∆D KL (t + 2) < 0, even though ∆D KL (t + 1) = 0, except in the stationary solution.