# How to demonstrate Benford's Law using Mathematica?

I tried to use Mathematica to demonstrate Benford's law. My problem may be with my understanding of the law. The code below shows how many times each digit from 1 through 9 is the leading digit in all integers from 2 to 10^7.

firstDigits=First[IntegerDigits[#]]&/@Range[2,10^7];
Counts[firstDigits]
(* <| 2->1111111,3->1111111,4->1111111,5->1111111,6->1111111,
7->1111111,8->1111111,9->1111111,1->1111111 |> *)


Above we see that 1, 2, 3, 4, 5, 6, 7, 8, 9 are each the leading digit exactly the same number of times. However, Benford's law says approximately 30.1% of the time the leading digit is 1. How can we use Mathematica to illustrate Benford's Law?

• Benford's law is not about just any old random collection of numbers (including the collection of all numbers, for which clearly each numbers shows up the same number of times as the leading digit). Benford's Law is about collections of numbers used for something, like the collection of telephone numbers, or house numbers, etc. If you read the page linked, it talks about Benford's law holding for data sets. Commented Aug 25, 2022 at 20:28
• AFAIR, scale invariance is an assumption for Benford's Law. E.g., firstDigits = First@First[RealDigits[#]] & /@ Exp@Range[-100., 100.]; Counts[firstDigits]/201. // KeySort Commented Aug 25, 2022 at 20:39
• You get a better approximation the closer you get to covering the positive reals at a fixed finite precision. RealDigits is slow, so here's a faster approach: firstDigits = 10^(# - UnitStep[#] + 1 &@ FractionalPart@Log10@Exp@Subdivide[-5000.24, 5000.24, 2^20]) // IntegerPart; Counts[firstDigits]/N@Length@firstDigits // KeySort Commented Aug 25, 2022 at 21:09
• There are examples on the Wolfram Demonstrations site. Commented Aug 25, 2022 at 21:20
• reference.wolfram.com/language/ref/BenfordDistribution.html Commented Aug 26, 2022 at 3:22

Here's an example using the population of all U.S. cities:

(* first get list of U.S. cities *)
cData = CityData[{All, "UnitedStates"}];
totalCities=Length@cData
(*32734 *)

(* retrieve population data from the city list *)
pVal = (QuantityMagnitude@CityData[#, "Population"]) & /@ cData;
firstDigits = First[IntegerDigits[#]] & /@ pVal;
(* tally the first digits in the list *)
tVals = Tally[firstDigits]
(* check the percentage of first digits that are 1 *)
oneDigits = Position[tVals, x_List /; x[[1]] == 1][[1, 1]]
oneCount = tVals[[oneDigits, 2]]
(* compute percent of  population that begins with digit 1*)
oneCount/totalCities // N
(* compare to the expected value *)
digitProb[n_] := Log[1 + 1/n]/Log[10];
digitProb[1] // N

(*0.302804*)

(* 0.30103*)