

 This page uses content from Wikipedia. The original article was at Benford's_law. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)


Benford's law, also called the first-digit law, refers to the frequency distribution of digits in many (but not all) real-life sources of data.

In this distribution, the number 1 occurs as the first digit about 30% of the time, while larger numbers occur in that position less frequently: 9 as the first digit less than 5% of the time. This distribution of first digits is the same as the widths of gridlines on a logarithmic scale.

Benford's law also concerns the expected distribution for digits beyond the first, which approach a uniform distribution.

This result has been found to apply to a wide variety of data sets, including electricity bills, street addresses, stock prices, population numbers, death rates, lengths of rivers, physical and mathematical constants, and processes described by power laws (which are very common in nature). It tends to be most accurate when values are distributed across multiple orders of magnitude.

A set of numbers is said to satisfy Benford's law if the leading digit 



d


{\displaystyle d}

  (



d
∈
{
1
,
…
,
9
}


{\displaystyle d\in \{1,\ldots ,9\}}

) occurs with probability

 



P
(
d
)
=

log

10


⁡
(
d
+
1
)
−

log

10


⁡
(
d
)
=

log

10


⁡

(
1
+


1
d


)



{\displaystyle P(d)=\log _{10}(d+1)-\log _{10}(d)=\log _{10}\left(1+{\frac {1}{d}}\right)}


For this task, write (a) routine(s) to calculate the distribution of first significant (non-zero) digits in a collection of numbers, then display the actual vs. expected distribution in the way most convenient for your language (table / graph / histogram / whatever).

Use the first 1000 numbers from the Fibonacci sequence as your data set. No need to show how the Fibonacci numbers are obtained.

You can generate them or load them from a file; whichever is easiest.

Display your actual vs expected distribution.


For extra credit: Show the distribution for one other set of numbers from a page on Wikipedia. State which Wikipedia page it can be obtained from and what the set enumerates. Again, no need to display the actual list of numbers or the code to load them.



See also:
 numberphile.com.
 A starting page on Wolfram Mathworld is Benfords Law .


# Generate the first n Fibonacci numbers: 1, 1, ...
# Numerical accuracy is insufficient beyond about 1450.
def fibonacci(n):
  # input: [f(i-2), f(i-1), countdown]
  def fib: (.[0] + .[1]) as $sum
           | if .[2] <= 0 then empty
             elif .[2] == 1 then $sum
             else $sum, ([ .[1], $sum, .[2] - 1 ] | fib)
             end;
  [1, 0, n] | fib ;
 
# is_prime is tailored to work with jq 1.4
def is_prime:
  if . == 2 then true
  else 2 < . and . % 2 == 1 and
       . as $in
       | (($in + 1) | sqrt) as $m
       | (((($m - 1) / 2) | floor) + 1) as $max
       | reduce range(1; $max) as $i
           (true; if . then ($in % ((2 * $i) + 1)) > 0 else false end)
  end ;
 
# primes in [m,n)
def primes(m;n):
  range(m;n) | select(is_prime);
 
def runs:
  reduce .[] as $item
    ( [];
      if . == [] then [ [ $item, 1] ] 
      else  .[length-1] as $last
            | if $last[0] == $item
              then (.[0:length-1] + [ [$item, $last[1] + 1] ] )
              else . + [[$item, 1]]
              end
      end ) ;
 
# Inefficient but brief:
def histogram: sort | runs;
 
def benford_probability:
  tonumber
  | if . > 0 then ((1 + (1 /.)) | log) / (10|log) 
    else 0
    end ;
 
# benford takes a stream and produces an array of [ "d", observed, expected ]
def benford(stream):
  [stream | tostring | .[0:1] ] | histogram as $histogram
  | reduce ($histogram | .[] | .[0]) as $digit
      ([]; . + [$digit, ($digit|benford_probability)] ) 
  | map(select(type == "number")) as $probabilities
  | ([ $histogram | .[] | .[1] ] | add) as $total
  | reduce range(0; $histogram|length) as $i
      ([]; . + ([$histogram[$i] + [$total * $probabilities[$i]] ] ) ) ;
 
# given an array of [value, observed, expected] values,
# produce the χ² statistic
def chiSquared:
  reduce .[] as $triple
    (0;
     if $triple[2] == 0 then .
     else . + ($triple[1] as $o | $triple[2] as $e | ($o - $e) | (.*.)/$e)
     end) ;
 
# truncate n places after the decimal point; 
# return a string since it can readily be converted back to a number
def precision(n):
  tostring as $s | $s | index(".")
  | if . then $s[0:.+n+1] else $s end ;
 
# Right-justify but do not truncate
def rjustify(n):
  length as $length | if n <= $length then . else " " * (n-$length) + . end;
 
# Attempt to align decimals so integer part is in a field of width n
def align(n):
  index(".") as $ix
  | if n < $ix then .
    elif $ix then (.[0:$ix]|rjustify(n)) +.[$ix:]
    else rjustify(n)
    end ;
 
# given an array of [value, observed, expected] values, 
# produce rows of the form: value observed expected
def print_rows(prec):
  .[] | map( precision(prec)|align(5) + "  ") | add ;
 
def report(heading; stream):
    benford(stream) as $array
    | heading,
      " Digit Observed Expected",
      ( $array | print_rows(2) ),
      "",
      " χ² = \( $array | chiSquared | precision(4))",
      ""
;
 
def task:
  report("First 100 fibonacci numbers:"; fibonacci( 100) ),
  report("First 1000 fibonacci numbers:"; fibonacci(1000) ),
  report("Primes less than 1000:"; primes(2;1000)),
  report("Primes between 1000 and 10000:"; primes(1000;10000)),
  report("Primes less than 100000:"; primes(2;100000)) 
;
 
task