Benford's law: Difference between revisions
Add ABC
No edit summary |
Not a robot (talk | contribs) (Add ABC) |
||
| (114 intermediate revisions by 47 users not shown) | |||
Line 11:
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
::::
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).
Line 19:
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 [[Fibonacci sequence|generate]] them or load them [http://www.
Display your actual vs expected distribution.
Line 31:
* A starting page on Wolfram Mathworld is {{Wolfram|Benfords|Law}}.
<br><br>
=={{header|8th}}==
<syntaxhighlight lang="8th">
: n:log10e ` 1 10 ln / ` ;
with: n
: n:log10 \ n -- n
ln log10e * ;
: benford \ x -- x
1 swap / 1+ log10 ;
: fibs \ xt n
swap >r
0.0 1.0 rot
( dup r@ w:exec tuck + ) swap times
2drop rdrop ;
var counts
: init
a:new ( 0 a:push ) 9 times counts ! ;
: leading \ n -- n
"%g" s:strfmt
0 s:@ '0 - nip ;
: bump-digit \ n --
1 swap
counts @ swap 1- ' + a:op! drop ;
: count-fibs \ --
( leading bump-digit ) 1000 fibs ;
: adjust \ --
counts @ ( 0.001 * ) a:map counts ! ;
: spaces \ n --
' space swap times ;
: .fw \ s n --
>r s:len r> rot . swap - spaces ;
: .header \ --
"Digit" 8 .fw "Expected" 10 .fw "Actual" 10 .fw cr ;
: .digit \ n --
>s 8 .fw ;
: .actual \ n --
"%.3f" s:strfmt 10 .fw ;
: .expected \ n --
"%.4f" s:strfmt 10 .fw ;
: report \ --
.header
counts @
( swap 1+ dup benford swap
.digit .expected .actual cr )
a:each drop ;
: benford-test
init count-fibs adjust report ;
;with
benford-test
bye
</syntaxhighlight>
{{out}}
<pre>
Digit Expected Actual
1 0.3010 0.301
2 0.1761 0.177
3 0.1249 0.125
4 0.0969 0.096
5 0.0792 0.080
6 0.0669 0.067
7 0.0580 0.056
8 0.0512 0.053
9 0.0458 0.045
</pre>
=={{header|11l}}==
{{trans|D}}
<syntaxhighlight lang="11l">F get_fibs()
V a = 1.0
V b = 1.0
[Float] r
L 1000
r [+]= a
(a, b) = (b, a + b)
R r
F benford(seq)
V freqs = [(0.0, 0.0)] * 9
V seq_len = 0
L(d) seq
I d != 0
freqs[String(d)[0].code - ‘1’.code][1]++
seq_len++
L(&f) freqs
f = (log10(1.0 + 1.0 / (L.index + 1)), f[1] / seq_len)
R freqs
print(‘#9 #9 #9’.format(‘Actual’, ‘Expected’, ‘Deviation’))
L(p) benford(get_fibs())
print(‘#.: #2.2% | #2.2% | #.4%’.format(L.index + 1, p[1] * 100, p[0] * 100, abs(p[1] - p[0]) * 100))</syntaxhighlight>
{{out}}
<pre>
Actual Expected Deviation
1: 30.10% | 30.10% | 0.0030%
2: 17.70% | 17.61% | 0.0909%
3: 12.50% | 12.49% | 0.0061%
4: 9.60% | 9.69% | 0.0910%
5: 8.00% | 7.92% | 0.0819%
6: 6.70% | 6.69% | 0.0053%
7: 5.60% | 5.80% | 0.1992%
8: 5.30% | 5.12% | 0.1847%
9: 4.50% | 4.58% | 0.0757%
</pre>
=={{header|ABC}}==
<syntaxhighlight lang="abc">HOW TO RETURN fibonacci.numbers n:
PUT 1, 1 IN a, b
PUT {} IN fibo
FOR i IN {1..n}:
INSERT a IN fibo
PUT b, a+b IN a, b
RETURN fibo
HOW TO RETURN digit.distribution nums:
PUT {} IN digits
FOR i IN {1..9}: PUT i IN digits["`i`"]
PUT {} IN dist
FOR i IN {1..9}: PUT 0 IN dist[i]
FOR n IN nums:
PUT digits["`n`"|1] IN digit
PUT dist[digit] + 1 IN dist[digit]
FOR i IN {1..9}:
PUT dist[i] / #nums IN dist[i]
RETURN dist
PUT digit.distribution fibonacci.numbers 1000 IN observations
WRITE "Digit"<<6, "Expected">>10, "Observed">>10/
FOR d IN {1..9}:
WRITE d<<6, ((10 log (1 + 1/d))>>10)|10, observations[d]>>10/</syntaxhighlight>
{{out}}
<pre>Digit Expected Observed
1 0.30102999 0.301
2 0.17609125 0.177
3 0.12493873 0.125
4 0.09691001 0.096
5 0.07918124 0.08
6 0.06694678 0.067
7 0.05799194 0.056
8 0.05115252 0.053
9 0.04575749 0.045</pre>
=={{header|Ada}}==
Line 36 ⟶ 198:
The program reads the Fibonacci-Numbers from the standard input. Each input line is supposed to hold N, followed by Fib(N).
<
procedure Benford is
Line 84 ⟶ 246:
Ada.Text_IO.New_Line;
end loop;
end Benford;</
{{out}}
Line 104 ⟶ 266:
Input is the list of primes below 100,000 from [http://www.mathsisfun.com/numbers/prime-number-lists.html]. Since each line in that file holds prime and only a prime, but no ongoing counter, we must slightly modify the program by commenting out a single line:
<
We can also edit out the declaration of the variable "Counter" ...or live with a compiler warning about never reading or assigning that variable.
Line 125 ⟶ 287:
=={{header|Aime}}==
<
sum(text a, text b)
{
Line 131 ⟶ 293:
integer e, f, n, r;
e =
f =
r = 0;
Line 163 ⟶ 325:
}
}
Line 172 ⟶ 334:
text a, b, w;
a = "0";
Line 182 ⟶ 344:
b = w;
c = w[0] - '0';
i += 1;
}
}
Line 198 ⟶ 360:
n = 1000;
fibs(f, n);
Line 212 ⟶ 370:
while (i < 9) {
i += 1;
o_form("%8d/p3d3w16//p3d3w16/\n", i, 100 * log10(1 + 1r / i), f[i] * m);
}
}</
{{out}}
<pre> expected found
Line 231 ⟶ 386:
8 5.115 5.300
9 4.575 4.5 </pre>
=={{header|ALGOL 68}}==
{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}}
Uses Algol 68G's LONG LONG INT which has programmer specifiable precision.
<syntaxhighlight lang="algol68">BEGIN
# set the number of digits for LONG LONG INT values #
PR precision 256 PR
# returns the probability of the first digit of each non-zero number in s #
PROC digit probability = ( []LONG LONG INT s )[]REAL:
BEGIN
[ 1 : 9 ]REAL result;
# count the number of times each digit is the first #
[ 1 : 9 ]INT count := ( 0, 0, 0, 0, 0, 0, 0, 0, 0 );
FOR i FROM LWB s TO UPB s DO
LONG LONG INT v := ABS s[ i ];
IF v /= 0 THEN
WHILE v > 9 DO v OVERAB 10 OD;
count[ SHORTEN SHORTEN v ] +:= 1
FI
OD;
# calculate the probability of each digit #
INT number of elements = ( UPB s + 1 ) - LWB s;
FOR i TO 9 DO
result[ i ] := IF number of elements = 0 THEN 0 ELSE count[ i ] / number of elements FI
OD;
result
END # digit probability # ;
# outputs the digit probabilities of some numbers and those expected by Benford's law #
PROC compare to benford = ( []REAL actual )VOID:
FOR i TO 9 DO
print( ( "Benford: ", fixed( log( 1 + ( 1 / i ) ), -7, 3 )
, " actual: ", fixed( actual[ i ], -7, 3 )
, newline
)
)
OD # compare to benford # ;
# generate 1000 fibonacci numbers #
[ 0 : 1000 ]LONG LONG INT fn;
fn[ 0 ] := 0;
fn[ 1 ] := 1;
FOR i FROM 2 TO UPB fn DO fn[ i ] := fn[ i - 1 ] + fn[ i - 2 ] OD;
# get the probabilities of each first digit of the fibonacci numbers and #
# compare to the probabilities expected by Benford's law #
compare to benford( digit probability( fn ) )
END</syntaxhighlight>
{{out}}
<pre>
Benford: 0.301 actual: 0.301
Benford: 0.176 actual: 0.177
Benford: 0.125 actual: 0.125
Benford: 0.097 actual: 0.096
Benford: 0.079 actual: 0.080
Benford: 0.067 actual: 0.067
Benford: 0.058 actual: 0.056
Benford: 0.051 actual: 0.053
Benford: 0.046 actual: 0.045
</pre>
=={{header|APL}}==
{{works with|Dyalog APL}}
<syntaxhighlight lang="apl">task←{
benf ← ≢÷⍨(⍳9)(+/∘.=)(⍎⊃∘⍕)¨
fibs ← (⊢,(+/¯2↑⊢))⍣998⊢1 1
exp ← 10⍟1+÷⍳9
obs ← benf fibs
⎕←'Expected Actual'⍪5⍕exp,[1.5]obs
}</syntaxhighlight>
{{out}}
<pre>Expected Actual
0.30103 0.30100
0.17609 0.17700
0.12494 0.12500
0.09691 0.09600
0.07918 0.08000
0.06695 0.06700
0.05799 0.05600
0.05115 0.05300
0.04576 0.04500</pre>
=={{header|Arturo}}==
<syntaxhighlight lang="arturo">fib: [a b]: <= [0 1]
do.times:998 -> [a b]: @[b 'fib ++ <= a+b]
leading: fib | map 'x -> first ~"|x|"
| tally
print "digit actual expected"
loop 1..9 'd ->
print [
pad.right ~"|d|" 8
pad.right ~"|leading\[d] // 1000|" 9
log 1 + 1//d 10
]</syntaxhighlight>
{{out}}
<pre>digit actual expected
1 0.301 0.3010299956639811
2 0.177 0.1760912590556812
3 0.125 0.1249387366082999
4 0.095 0.09691001300805641
5 0.08 0.0791812460476248
6 0.067 0.06694678963061322
7 0.056 0.05799194697768673
8 0.053 0.05115252244738128
9 0.045 0.04575749056067514</pre>
=={{header|AutoHotkey}}==
{{works with|AutoHotkey_L}}(AutoHotkey1.1+)
<
fib := NStepSequence(1, 1, 2, 1000)
Out := "Digit`tExpected`tObserved`tDeviation`n"
Line 276 ⟶ 540:
}
return, (Carry ? Carry : "") . Result
}</
'''Output:'''
<pre>Digit Expected Observed Deviation
Line 290 ⟶ 554:
=={{header|AWK}}==
<syntaxhighlight lang="awk">
# syntax: GAWK -f BENFORDS_LAW.AWK
BEGIN {
Line 318 ⟶ 582:
function abs(x) { if (x >= 0) { return x } else { return -x } }
function log10(x) { return log(x)/log(10) }
</syntaxhighlight>
{{out}}
<pre>
Line 332 ⟶ 596:
9 4.5757 4.5000 0.0757
</pre>
=={{header|BCPL}}==
BCPL doesn't do floating point well, so I use integer math to compute the most significant digits of the Fibonacci sequence and use a table that has the values of log10(d + 1/d)
<syntaxhighlight lang="bcpl">
GET "libhdr"
MANIFEST {
fib_sz = 1_000_000_000 // keep 9 significant digits for computing F(n)
}
LET msd(n) = VALOF {
UNTIL n < 10 DO
n /:= 10
RESULTIS n
}
LET fibonacci(n, tally) BE {
LET a, b, c, e = 0, 1, ?, 0
FOR i = 1 TO n {
TEST e = 0
THEN tally[msd(b)] +:= 1
ELSE tally[b / (fib_sz / 10)] +:= 1
c := a + b
IF c > fib_sz {
a := a / 10 - (a MOD 10 >= 5) // subtract, since condition evalutes to
b := b / 10 - (b MOD 10 >= 5) // eithr 0 or -1.
c := a + b
e +:= 1 // keep track of exponent, just 'cuz
}
a, b := b, c
}
}
LET start() = VALOF {
// expected value of benford: log10(d + 1/d)
LET expected = TABLE 0, 301, 176, 125, 97, 79, 67, 58, 51, 46
LET actual = VEC 9
FOR i = 0 TO 9 DO actual!i := 0
fibonacci(1000, actual)
writes("*nLeading digital distribution of the first 1,000 Fibonacci numbers*n")
writes("Digit*tActual*tExpected*n")
FOR i = 1 TO 9
writef("%i *t%0.3d *t%0.3d *n", i, actual!i, expected!i)
RESULTIS 0
}
</syntaxhighlight>
{{Out}}
<pre>
BCPL 32-bit Cintcode System (13 Jan 2020)
0.000>
Leading digital distribution of the first 1,000 Fibonacci numbers
Digit Actual Expected
1 0.301 0.301
2 0.177 0.176
3 0.125 0.125
4 0.096 0.097
5 0.080 0.079
6 0.067 0.067
7 0.056 0.058
8 0.053 0.051
9 0.045 0.046
</pre>
=={{header|BASIC256}}==
<syntaxhighlight lang="vb">n = 1000
dim actual(n) fill 0
for nr = 1 to n
num$ = string(fibonacci(nr))
j = int(left(num$,1))
actual[j] += 1
next
print "First 1000 Fibonacci numbers"
print "Digit ", "Actual freq ", "Expected freq"
for i = 1 to 9
freq = frequency(i)*100
print " "; ljust(i,4), rjust(actual[i]/10,5), rjust(freq,5)
next
end
function frequency(n)
return (log10(n+1) - log10(n))
end function
function fibonacci(f)
f = int(f)
a = 0 : b = 1 : c = 0 : n = 0
while n < f
a = b
b = c
c = a + b
n += 1
end while
return c
end function</syntaxhighlight>
{{Out}}
<pre>First 1000 Fibonacci numbers
Digit Actual freq Expected freq
1 30.1 30.1029995664
2 17.7 17.6091259056
3 12.5 12.4938736608
4 9.6 9.69100130081
5 8.0 7.91812460476
6 6.7 6.694679
7 5.6 5.79919469777
8 5.3 5.11525224474
9 4.5 4.57574905607</pre>
=={{header|C}}==
<
#include <stdlib.h>
#include <math.h>
Line 397 ⟶ 770:
return EXIT_SUCCESS;
}</
{{Out}}
Line 412 ⟶ 785:
8 0.053 0.051
9 0.045 0.046</pre>
=={{header|C++}}==
<
//if used prepackaged you can compile writing "g++ -std=c++11 -lcln yourprogram.cpp -o yourprogram"
#include <cln/integer.h>
Line 476 ⟶ 848:
return 0 ;
}
</syntaxhighlight>
{{out}}
<pre> found expected
Line 490 ⟶ 862:
</pre>
=={{header|Chipmunk Basic}}==
{{trans|Yabasic}}
{{works with|Chipmunk Basic|3.6.4}}
<syntaxhighlight lang="vbnet">100 cls
110 n = 1000
120 dim actual(n)
130 for nr = 1 to n
140 num$ = str$(fibonacci(nr))
150 j = val(left$(num$,1))
160 actual(j) = actual(j)+1
170 next
180 print "First 1000 Fibonacci numbers"
190 print "Digit Actual freq Expected freq"
200 for i = 1 to 9
210 freq = frequency(i)*100
220 print format$(i,"###");
230 print using " ##.###";actual(i)/10;
240 print using " ##.###";freq
250 next
260 end
270 sub frequency(n)
280 frequency = (log10(n+1)-log10(n))
290 end sub
300 sub log10(n)
310 log10 = log(n)/log(10)
320 end sub
330 sub fibonacci(n)
335 rem https://rosettacode.org/wiki/Fibonacci_sequence#Chipmunk_Basic
340 n1 = 0
350 n2 = 1
360 for k = 1 to abs(n)
370 sum = n1+n2
380 n1 = n2
390 n2 = sum
400 next k
410 if n < 0 then
420 fibonacci = n1*((-1)^((-n)+1))
430 else
440 fibonacci = n1
450 endif
460 end sub</syntaxhighlight>
=={{header|Clojure}}==
<syntaxhighlight lang="lisp">(ns example
(:gen-class))
(defn abs [x]
(if (> x 0)
x
(- x)))
(defn calc-benford-stats [digits]
" Frequencies of digits in data "
(let [y (frequencies digits)
tot (reduce + (vals y))]
[y tot]))
(defn show-benford-stats [v]
" Prints in percent the actual, Benford expected, and difference"
(let [fd (map (comp first str) v)] ; first digit of each record
(doseq [q (range 1 10)
:let [[y tot] (calc-benford-stats fd)
d (first (str q)) ; reference digit
f (/ (get y d 0) tot 0.01) ; percent of occurence of digit
p (* (Math/log10 (/ (inc q) q)) 100) ; Benford expected percent
e (abs (- f p))]] ; error (difference)
(println (format "%3d %10.2f %10.2f %10.2f"
q
f
p
e)))))
; Generate fibonacci results
(def fib (lazy-cat [0N 1N] (map + fib (rest fib))))
;(def fib-digits (map (comp first str) (take 10000 fib)))
(def fib-digits (take 10000 fib))
(def header " found-% expected-% diff")
(println "Fibonacci Results")
(println header)
(show-benford-stats fib-digits)
;
; Universal Constants from Physics (using first column of data)
(println "Universal Constants from Physics")
(println header)
(let [
data-parser (fn [s]
(let [x (re-find #"\s{10}-?[0|/\.]*([1-9])" s)]
(if (not (nil? x)) ; Skips records without number
(second x)
x)))
input (slurp "http://physics.nist.gov/cuu/Constants/Table/allascii.txt")
y (for [line (line-seq (java.io.BufferedReader.
(java.io.StringReader. input)))]
(data-parser line))
z (filter identity y)]
(show-benford-stats z))
; Sunspots
(println "Sunspots average count per month since 1749")
(println header)
(let [
data-parser (fn [s]
(nth (re-find #"(.+?\s){3}([1-9])" s) 2))
; Sunspot data loaded from file (saved from ;https://solarscience.msfc.nasa.gov/greenwch/SN_m_tot_V2.0.txt")
; (note: attempting to load directly from url causes https Trust issues, so saved to file after loading to Browser)
input (slurp "SN_m_tot_V2.0.txt")
y (for [line (line-seq (java.io.BufferedReader.
(java.io.StringReader. input)))]
(data-parser line))]
(show-benford-stats y))
</syntaxhighlight>
{{Output}}
<pre>
Fibonacci Results
found-% expected-% diff
1 30.11 30.10 0.01
2 17.62 17.61 0.01
3 12.49 12.49 0.00
4 9.68 9.69 0.01
5 7.92 7.92 0.00
6 6.68 6.69 0.01
7 5.80 5.80 0.00
8 5.13 5.12 0.01
9 4.56 4.58 0.02
Universal Constants from Physics
found-% expected-% diff
1 34.34 30.10 4.23
2 18.67 17.61 1.07
3 9.04 12.49 3.46
4 8.43 9.69 1.26
5 8.43 7.92 0.52
6 7.23 6.69 0.53
7 3.31 5.80 2.49
8 5.12 5.12 0.01
9 5.42 4.58 0.85
Sunspots average count per month since 1749
found-% expected-% diff
1 37.44 30.10 7.34
2 16.28 17.61 1.33
3 7.16 12.49 5.34
4 6.88 9.69 2.81
5 6.35 7.92 1.57
6 6.04 6.69 0.66
7 7.25 5.80 1.45
8 5.57 5.12 0.46
9 5.76 4.58 1.18
</pre>
=={{header|CoffeeScript}}==
<
a = 1; b = 0
return () ->
Line 511 ⟶ 1,037:
console.log "Digit\tActual\tExpected"
for i in [1..9]
console.log i + "\t" + actual[i - 1].toFixed(3) + '\t' + expected[i - 1].toFixed(3)</
{{out}}
<pre>Leading digital distribution of the first 1,000 Fibonacci numbers
Line 524 ⟶ 1,050:
8 0.053 0.051
9 0.045 0.046</pre>
=={{header|Common Lisp}}==
<
"Return the frequency distribution of the most significant nonzero
digits in the given list of numbers. The first element of the list
Line 568 ⟶ 1,093:
(map 'list #'list '(1 2 3 4 5 6 7 8 9)
actual-distribution
expected-distribution))))</
<pre>; *fib1000* is a list containing the first 1000 numbers in the Fibonnaci sequence
Line 582 ⟶ 1,107:
8 0.053 0.051
9 0.045 0.046</pre>
=={{header|Crystal}}==
{{trans|Ruby}}
<syntaxhighlight lang="ruby">require "big"
EXPECTED = (1..9).map{ |d| Math.log10(1 + 1.0 / d) }
def fib(n)
a, b = 0.to_big_i, 1.to_big_i
(0...n).map { ret, a, b = a, b, a + b; ret }
end
# powers of 3 as a test sequence
def power_of_threes(n)
(0...n).map { |k| 3.to_big_i ** k }
end
def heads(s)
s.map { |a| a.to_s[0].to_i }
end
def show_dist(title, s)
s = heads(s)
c = Array.new(10, 0)
s.each{ |x| c[x] += 1 }
siz = s.size
res = (1..9).map{ |d| c[d] / siz }
puts "\n %s Benfords deviation" % title
res.zip(EXPECTED).each_with_index(1) do |(r, e), i|
puts "%2d: %5.1f%% %5.1f%% %5.1f%%" % [i, r*100, e*100, (r - e).abs*100]
end
end
def random(n)
(0...n).map { |i| rand(1..n) }
end
show_dist("fibbed", fib(1000))
show_dist("threes", power_of_threes(1000))
# just to show that not all kind-of-random sets behave like that
show_dist("random", random(10000))</syntaxhighlight>
{{out}}
<pre>
fibbed Benfords deviation
1: 30.1% 30.1% 0.0%
2: 17.7% 17.6% 0.1%
3: 12.5% 12.5% 0.0%
4: 9.5% 9.7% 0.2%
5: 8.0% 7.9% 0.1%
6: 6.7% 6.7% 0.0%
7: 5.6% 5.8% 0.2%
8: 5.3% 5.1% 0.2%
9: 4.5% 4.6% 0.1%
threes Benfords deviation
1: 30.0% 30.1% 0.1%
2: 17.7% 17.6% 0.1%
3: 12.3% 12.5% 0.2%
4: 9.8% 9.7% 0.1%
5: 7.9% 7.9% 0.0%
6: 6.6% 6.7% 0.1%
7: 5.9% 5.8% 0.1%
8: 5.2% 5.1% 0.1%
9: 4.6% 4.6% 0.0%
random Benfords deviation
1: 11.2% 30.1% 18.9%
2: 10.8% 17.6% 6.8%
3: 10.8% 12.5% 1.7%
4: 11.0% 9.7% 1.3%
5: 11.1% 7.9% 3.1%
6: 10.8% 6.7% 4.1%
7: 10.8% 5.8% 5.0%
8: 11.2% 5.1% 6.1%
9: 12.3% 4.6% 7.7%
</pre>
=={{header|D}}==
{{trans|Scala}}
<
double[2][9] benford(R)(R seq) if (isForwardRange!R && !isInfinite!R) {
Line 608 ⟶ 1,209:
writefln("%d: %5.2f%% | %5.2f%% | %5.4f%%",
i+1, p[1] * 100, p[0] * 100, abs(p[1] - p[0]) * 100);
}</
{{out}}
<pre> Actual Expected Deviation
Line 624 ⟶ 1,225:
===Alternative Version===
The output is the same.
<
std.algorithm, std.array;
Line 641 ⟶ 1,242:
f * 100.0 / N, expected[i] * 100,
abs((f / double(N)) - expected[i]) * 100);
}</
=={{header|Delphi}}==
See [[Benford's law#Pascal|Pascal]].
=={{header|EasyLang}}==
<syntaxhighlight lang=easylang>
func$ add a$ b$ .
for i to higher len a$ len b$
a = number substr a$ i 1
b = number substr b$ i 1
r = a + b + c
c = r div 10
r$ &= r mod 10
.
if c > 0
r$ &= c
.
return r$
.
#
len fibdist[] 9
proc mkfibdist . .
# generate 1000 fibonacci numbers as
# (reversed) strings, because 53 bit
# integers are too small
#
n = 1000
prev$ = 0
val$ = 1
fibdist[1] = 1
for i = 2 to n
h$ = add prev$ val$
prev$ = val$
val$ = h$
ind = number substr val$ len val$ 1
fibdist[ind] += 1
.
for i to len fibdist[]
fibdist[i] = fibdist[i] / n
.
.
mkfibdist
#
len benfdist[] 9
proc mkbenfdist . .
for i to 9
benfdist[i] = log10 (1 + 1.0 / i)
.
.
mkbenfdist
#
numfmt 3 0
print "Actual Expected"
for i to 9
print fibdist[i] & " " & benfdist[i]
.
</syntaxhighlight>
{{out}}
<pre>
Actual Expected
0.301 0.301
0.177 0.176
0.125 0.125
0.096 0.097
0.080 0.079
0.067 0.067
0.056 0.058
0.053 0.051
0.045 0.046
</pre>
=={{header|Elixir}}==
<
def distribution(n), do: :math.log10( 1 + (1 / n) )
Line 661 ⟶ 1,331:
end
Benfords_law.task</
{{out}}
Line 678 ⟶ 1,348:
=={{header|Erlang}}==
<syntaxhighlight lang="erlang">
-module( benfords_law ).
-export( [actual_distribution/1, distribution/1, task/0] ).
Line 704 ⟶ 1,374:
[Key | _] = erlang:integer_to_list( N ),
dict:update_counter( Key - 48, 1, Dict ).
</syntaxhighlight>
{{out}}
<pre>
Line 719 ⟶ 1,389:
9 0.045 0.04575749056067514
</pre>
=={{header|F sharp|F#}}==
For Fibonacci code, see https://rosettacode.org/wiki/Fibonacci_sequence#F.89
<syntaxhighlight lang="fsharp">open System
let fibonacci = Seq.unfold (fun (x, y) -> Some(x, (y, x + y))) (0I,1I)
let fibFirstNumbers nth =
fibonacci |> Seq.take nth |> Seq.map (fun n -> n.ToString().[0] |> string |> Int32.Parse)
let fibFirstNumbersFrequency nth =
let firstNumbers = fibFirstNumbers nth |> Seq.toList
let counts = firstNumbers |> List.countBy id |> List.sort |> List.filter (fun (k, _) -> k <> 0)
let total = firstNumbers |> List.length |> float
counts |> List.map (fun (_, v) -> float v/total)
let benfordLaw d = log10(1.0 + (1.0/float d))
let benfordLawFigures = [1..9] |> List.map benfordLaw
let run () =
printfn "Frequency of the first digit (1 through 9) in the Fibonacci sequence:"
fibFirstNumbersFrequency 1000 |> List.iter (fun f -> printf $"{f:N5} ")
printfn "\nBenford's law for 1 through 9:"
benfordLawFigures |> List.iter (fun f -> printf $"{f:N5} ")
</syntaxhighlight>
{{out}}
<pre>
Frequency of the first digit (1 through 9) in the Fibonacci sequence:
0.30100 0.17700 0.12500 0.09500 0.08000 0.06700 0.05600 0.05300 0.04500
Benford's law for 1 through 9:
0.30103 0.17609 0.12494 0.09691 0.07918 0.06695 0.05799 0.05115 0.04576
</pre>
=={{header|Factor}}==
<syntaxhighlight lang="factor">USING: assocs compiler.tree.propagation.call-effect formatting
kernel math math.functions math.statistics math.text.utils
sequences ;
IN: rosetta-code.benfords-law
: expected ( n -- x ) recip 1 + log10 ;
: next-fib ( vec -- vec' )
[ last2 ] keep [ + ] dip [ push ] keep ;
: data ( -- seq ) V{ 1 1 } clone 998 [ next-fib ] times ;
: 1st-digit ( n -- m ) 1 digit-groups last ;
: leading ( -- seq ) data [ 1st-digit ] map ;
: .header ( -- )
"Digit" "Expected" "Actual" "%-10s%-10s%-10s\n" printf ;
: digit-report ( digit digit-count -- digit expected actual )
dupd [ expected ] dip 1000 /f ;
: .digit-report ( digit digit-count -- )
digit-report "%-10d%-10.4f%-10.4f\n" printf ;
: main ( -- )
.header leading histogram [ .digit-report ] assoc-each ;
MAIN: main</syntaxhighlight>
{{out}}
<pre>
Digit Expected Actual
1 0.3010 0.3010
2 0.1761 0.1770
3 0.1249 0.1250
4 0.0969 0.0960
5 0.0792 0.0800
6 0.0669 0.0670
7 0.0580 0.0560
8 0.0512 0.0530
9 0.0458 0.0450
</pre>
=={{header|Forth}}==
<
: f2drop fdrop fdrop ;
Line 764 ⟶ 1,510:
set-precision ;
: compute-benford tally report ;</
{{Out}}
<pre>Gforth 0.7.2, Copyright (C) 1995-2008 Free Software Foundation, Inc.
Line 781 ⟶ 1,527:
8 0.053 0.0512
9 0.045 0.0458 ok</pre>
=={{header|Fortran}}==
FORTRAN 90. Compilation and output of this program using emacs compile command and a fairly obvious Makefile entry:
<
Compilation started at Sat May 18 01:13:00
Line 792 ⟶ 1,537:
0.300999999 0.177000001 0.125000000 9.60000008E-02 7.99999982E-02 6.70000017E-02 5.70000000E-02 5.29999994E-02 4.50000018E-02 LEADING FIBONACCI DIGIT
Compilation finished at Sat May 18 01:13:00</
<
! compute most significant digits, Fibonacci like.
implicit none
Line 862 ⟶ 1,607:
write(6,*) (benfordsLaw(i),i=1,9),'THE LAW'
write(6,*) (count(i)/1000.0 ,i=1,9),'LEADING FIBONACCI DIGIT'
end program benford</
=={{header|FreeBASIC}}==
{{libheader|GMP}}
<syntaxhighlight lang="freebasic">' version 27-10-2016
' compile with: fbc -s console
#Define max 1000 ' total number of Fibonacci numbers
#Define max_sieve 15485863 ' should give 1,000,000
#Include Once "gmp.bi" ' uses the GMP libary
Dim As ZString Ptr z_str
Dim As ULong n, d
ReDim As ULong digit(1 To 9)
Dim As Double expect, found
Dim As mpz_ptr fib1, fib2
fib1 = Allocate(Len(__mpz_struct)) : Mpz_init_set_ui(fib1, 0)
fib2 = Allocate(Len(__mpz_struct)) : Mpz_init_set_ui(fib2, 1)
digit(1) = 1 ' fib2
For n = 2 To max
Swap fib1, fib2 ' fib1 = 1, fib2 = 0
mpz_add(fib2, fib1, fib2) ' fib1 = 1, fib2 = 1 (fib1 + fib2)
z_str = mpz_get_str(0, 10, fib2)
d = Val(Left(*z_str, 1)) ' strip the 1 digit on the left off
digit(d) = digit(d) +1
Next
mpz_clear(fib1) : DeAllocate(fib1)
mpz_clear(fib2) : DeAllocate(fib2)
Print
Print "First 1000 Fibonacci numbers"
Print "nr: total found expected difference"
For d = 1 To 9
n = digit(d)
found = n / 10
expect = (Log(1 + 1 / d) / Log(10)) * 100
Print Using " ## ##### ###.## % ###.## % ##.### %"; _
d; n ; found; expect; expect - found
Next
ReDim digit(1 To 9)
ReDim As UByte sieve(max_sieve)
'For d = 4 To max_sieve Step 2
' sieve(d) = 1
'Next
Print : Print "start sieve"
For d = 3 To sqr(max_sieve)
If sieve(d) = 0 Then
For n = d * d To max_sieve Step d * 2
sieve(n) = 1
Next
End If
Next
digit(2) = 1 ' 2
Print "start collecting first digits"
For n = 3 To max_sieve Step 2
If sieve(n) = 0 Then
d = Val(Left(Trim(Str(n)), 1))
digit(d) = digit(d) +1
End If
Next
Dim As ulong total
For n = 1 To 9
total = total + digit(n)
Next
Print
Print "First";total; " primes"
Print "nr: total found expected difference"
For d = 1 To 9
n = digit(d)
found = n / total * 100
expect = (Log(1 + 1 / d) / Log(10)) * 100
Print Using " ## ######## ###.## % ###.## % ###.### %"; _
d; n ; found; expect; expect - found
Next
' empty keyboard buffer
While InKey <> "" : Wend
Print : Print "hit any key to end program"
Sleep
End</syntaxhighlight>
{{out}}
<pre>First 1000 Fibonacci numbers
nr: total found expected difference
1 301 30.10 % 30.10 % 0.003 %
2 177 17.70 % 17.61 % -0.091 %
3 125 12.50 % 12.49 % -0.006 %
4 96 9.60 % 9.69 % 0.091 %
5 80 8.00 % 7.92 % -0.082 %
6 67 6.70 % 6.69 % -0.005 %
7 56 5.60 % 5.80 % 0.199 %
8 53 5.30 % 5.12 % -0.185 %
9 45 4.50 % 4.58 % 0.076 %
start sieve
start collecting first digits
First1000000 primes
nr: total found expected difference
1 415441 41.54 % 30.10 % -11.441 %
2 77025 7.70 % 17.61 % 9.907 %
3 75290 7.53 % 12.49 % 4.965 %
4 74114 7.41 % 9.69 % 2.280 %
5 72951 7.30 % 7.92 % 0.623 %
6 72257 7.23 % 6.69 % -0.531 %
7 71564 7.16 % 5.80 % -1.357 %
8 71038 7.10 % 5.12 % -1.989 %
9 70320 7.03 % 4.58 % -2.456 %</pre>
=={{header|Fōrmulæ}}==
{{FormulaeEntry|page=https://formulae.org/?script=examples/Benford%27s_law}}
'''Solution:'''
The following function calculates the distribution of the first digit of a list of integer, positive numbers:
[[File:Fōrmulæ - Benford's law 01.png]]
'''Example:'''
[[File:Fōrmulæ - Benford's law 02.png]]
[[File:Fōrmulæ - Benford's law 03.png]]
'''Benford distribution.''' Benford distribution is, by definition (using a precision of 10 digits):
[[File:Fōrmulæ - Benford's law 04.png]]
[[File:Fōrmulæ - Benford's law 05.png]]
'''Comparison chart.''' The following function creates a chart, in order to compare the distribution of the first digis in a givel list or numbers, and the Benford distribution:
[[File:Fōrmulæ - Benford's law 06.png]]
'''Testing.''' Testing with the previous list of numbers:
[[File:Fōrmulæ - Benford's law 07.png]]
[[File:Fōrmulæ - Benford's law 08.png]]
'''Case 1.''' Use the first 1,000 numbers from the Fibonacci sequence as your data set
[[File:Fōrmulæ - Benford's law 09.png]]
[[File:Fōrmulæ - Benford's law 10.png]]
'''Case 2.''' The sequence of the fisrt 1,000 natural numbers does not follow the Benford distribution:
[[File:Fōrmulæ - Benford's law 11.png]]
[[File:Fōrmulæ - Benford's law 12.png]]
'''Case 3.''' The following example is for the list of the fist 1,000 factorial numbers.
[[File:Fōrmulæ - Benford's law 13.png]]
[[File:Fōrmulæ - Benford's law 14.png]]
=={{header|FutureBasic}}==
<syntaxhighlight lang="futurebasic">
Short t, i, j, k, m
Double a(9), z
Double phi, psi
CFStringRef s
print @"Benford:"
for i = 1 to 9
a(i) = log10( 1 + 1 / i )
print fn StringWithFormat( @"%.3f ", a(i) ),
next
// Fibonacci according to DeMoivre and Binet
for t = 1 to 9 : a(t) = 0 : next // Clean array
phi = ( 1 + sqr(5) ) / 2
psi = ( 1 - sqr(5) ) / 2
for i = 1 to 1000
z = ( phi^i - psi^i ) / sqr( 5 )
s = fn StringWithFormat( @"%e", z) // Get first digit
t = fn StringIntegerValue( left( s, 1 ) )
a(t) = a(t) + 1
next
print @"\n\nFibonacci:"
for i = 1 to 9
print fn StringWithFormat( @"%.3f ", a(i) / 1000 ),
next
// Multiplication tables
for t = 1 to 9 : a(t) = 0 : next // Clean array
for i = 1 to 10
for j = 1 to 10
for k = 1 to 10
for m = 1 to 10
z = i * j * k * m
s = fn StringWithFormat( @"%e", z )
t = fn StringIntegerValue( left( s, 1 ) )
a(t) = a(t) + 1
next
next
next
next
print @"\n\nMultiplication:"
for i = 1 to 9
print fn StringWithFormat( @"%.3f ", a(i) / 1e4 ),
next
// Factorials according to DeMoivre and Stirling
for t = 1 to 9 : a(t) = 0 : next // Clean array
for i = 10 to 110
z = sqr(2 * pi * i ) * (i / exp(1) )^i
s = fn StringWithFormat( @"%e", z )
t = fn StringIntegerValue( left( s, 1 ) )
a(t) = a(t) + 1
next
print @"\n\nFactorials:"
for i = 1 to 9
print fn StringWithFormat( @"%.2f ", a(i) / 100 ),
next
handleevents
}</syntaxhighlight>
{{Out}}
<pre>
Benford:
0.301 0.176 0.125 0.097 0.079 0.067 0.058 0.051 0.046
Fibonacci:
0.301 0.177 0.125 0.096 0.080 0.067 0.056 0.053 0.045
Multiplication:
0.302 0.181 0.124 0.103 0.071 0.061 0.055 0.055 0.047
Factorials:
0.35 0.16 0.12 0.06 0.06 0.06 0.02 0.10 0.08
</pre>
=={{header|Go}}==
<
import (
Line 895 ⟶ 1,891:
math.Log10(1+1/float64(i+1)))
}
}</
{{Out}}
<pre>
Line 909 ⟶ 1,905:
8 0.053 0.051
9 0.045 0.046
</pre>
=={{header|Groovy}}==
'''Solution:'''<br>
Uses [[Fibonacci_sequence#Analytic_8|Fibonacci sequence analytic formula]]
{{trans|Java}}
<syntaxhighlight lang="groovy">def tallyFirstDigits = { size, generator ->
def population = (0..<size).collect { generator(it) }
def firstDigits = [0]*10
population.each { number ->
firstDigits[(number as String)[0] as int] ++
}
firstDigits
}</syntaxhighlight>
'''Test:'''
<syntaxhighlight lang="groovy">def digitCounts = tallyFirstDigits(1000, aFib)
println "d actual predicted"
(1..<10).each {
printf ("%d %10.6f %10.6f\n", it, digitCounts[it]/1000, Math.log10(1.0 + 1.0/it))
}</syntaxhighlight>
'''Output:'''
<pre>d actual predicted
1 0.301000 0.301030
2 0.177000 0.176091
3 0.125000 0.124939
4 0.095000 0.096910
5 0.080000 0.079181
6 0.067000 0.066947
7 0.056000 0.057992
8 0.053000 0.051153
9 0.045000 0.045757</pre>
=={{header|GW-BASIC}}==
{{improve|GW-BASIC|Does not show actual vs expected values for Fibonaccis as task specifies.}}
<syntaxhighlight lang="GW-BASIC">
10 DEF FNBENFORD(N)=LOG(1+1/N)/LOG(10)
20 CLS
30 PRINT "One digit Benford's Law"
40 FOR i = 1 TO 9
50 PRINT i,:PRINT USING "##.######";FNBENFORD(i)
60 NEXT i
70 END
</syntaxhighlight>
{{out}}
<pre>
1 0.301030
2 0.176091
3 0.124939
4 0.096910
5 0.079181
6 0.066947
7 0.057992
8 0.051153
9 0.045758
</pre>
=={{header|Haskell}}==
<
import Data.Char (digitToInt)
Line 918 ⟶ 1,969:
fstdigit = digitToInt . head . show
n = 1000 :: Int
fibs = 1 : 1 : zipWith (+) fibs (tail fibs)
fibdata = map fstdigit $ take n fibs
freqs = M.fromListWith (+) $ zip fibdata (repeat 1)
tab ::
tab =
[ ( d
, fromIntegral (M.findWithDefault
, logBase 10.0 $ 1 + 1 / fromIntegral d)
| d <- [1 .. 9] ]
main = print tab</
{{out}}
<pre>[(1,0.301,0.301029995663981),
Line 940 ⟶ 1,996:
(9,0.045,0.0457574905606751)]
</pre>
=={{header|Icon}} and {{header|Unicon}}==
The following solution works in both languages.
<
procedure main()
Line 978 ⟶ 2,033:
if n%2 = 1 then return [c+d, d]
else return [d, c]
end</
Sample run:
Line 995 ⟶ 2,050:
->
</pre>
=={{header|J}}==
'''Solution'''
<syntaxhighlight lang="j">benford=: 10&^.@(1 + %) NB. expected frequency of first digit y
Digits=: '123456789'
firstSigDigits=: {.@(-. -.&Digits)@":"0 NB. extract first significant digit from numbers
freq=: (] % +/)@:<:@(#/.~)@, NB. calc frequency of values (x) in y</syntaxhighlight>
'''Required Example'''
<syntaxhighlight lang="j"> First1000Fib=: (, +/@:(_2&{.)) ^: (1000-#) 1 1
NB. Expected vs Actual frequencies for Digits 1-9
Digits ((] ,. benford)@"."0@[ ,. (freq firstSigDigits)) First1000Fib
1 0.30103 0.301
2 0.176091 0.177
3 0.124939 0.125
4 0.09691 0.096
5 0.0791812 0.08
6 0.0669468 0.067
7 0.0579919 0.056
8 0.0511525 0.053
9 0.0457575 0.045</syntaxhighlight>
'''Alternatively'''
We show the correlation coefficient of Benford's law with the leading digits of the first 1000 Fibonacci numbers is almost unity.
<
benford =: log10@:(1+%)
assert '0.30 0.18 0.12 0.10 0.08 0.07 0.06 0.05 0.05' -: 5j2 ": benford >: i. 9
Line 1,033 ⟶ 2,109:
assert '0.9999' -: 6j4 ": (normalize TALLY_BY_KEY) r benford >: i.9
assert '0.9999' -: 6j4 ": TALLY_BY_KEY r benford >: i.9 NB. Of course we don't need normalization</
=={{header|Java}}==
<
import java.util.Locale;
public class
private static
}
return fib;
}
public static void main(String[] args) {
BigInteger[] numbers = generateFibonacci(1000);
int[] firstDigits = new int[10];
firstDigits[Integer.valueOf(number.toString().substring(0, 1))]++;
}
for (int i = 1; i < firstDigits.length; i++) {
System.out.printf(Locale.ROOT, "%d %10.6f %10.6f%n",
}
}</syntaxhighlight>
The output is:
<pre>1
2 0.177000 0.176091
3 0.125000 0.124939
4 0.096000 0.096910
5 0.080000 0.079181
6 0.067000 0.066947
7 0.056000 0.057992
8 0.053000 0.051153
9 0.045000 0.045757</pre>
To use other number sequences, implement a suitable <tt>NumberGenerator</tt>, construct a <tt>Benford</tt> instance with it and print it.
=={{header|JavaScript}}==
<syntaxhighlight lang="javascript">const fibseries = n => [...Array(n)]
.reduce(
(fib, _, i) => i < 2 ? (
fib
) : fib.concat(fib[i - 1] + fib[i - 2]),
[1, 1]
);
const benford = array => [1, 2, 3, 4, 5, 6, 7, 8, 9]
.map(val => [val, array
.reduce(
(sum, item) => sum + (
`${item}` [0] === `${val}`
),
0
) / array.length, Math.log10(1 + 1 / val)
]);
console.log(benford(fibseries(1000)))</syntaxhighlight>
{{output}}
<pre>0: (3) [1, 0.301, 0.3010299956639812]
1: (3) [2, 0.177, 0.17609125905568124]
2: (3) [3, 0.125, 0.12493873660829992]
3: (3) [4, 0.096, 0.09691001300805642]
4: (3) [5, 0.08, 0.07918124604762482]
5: (3) [6, 0.067, 0.06694678963061322]
6: (3) [7, 0.056, 0.05799194697768673]
7: (3) [8, 0.053, 0.05115252244738129]
8: (3) [9, 0.045, 0.04575749056067514]</pre>
=={{header|jq}}==
{{works with|jq|1.4}}
This implementation shows the observed and expected number of occurrences together with the χ² statistic.
For the sake of being self-contained, the following includes a generator for Fibonacci numbers, and a prime number generator that is inefficient but brief and can generate numbers within an arbitrary range.<
# Numerical accuracy is insufficient beyond about 1450.
def fibonacci(n):
Line 1,199 ⟶ 2,292:
;
task</
{{out}}
<pre>First 100 fibonacci numbers:
Line 1,271 ⟶ 2,364:
χ² = 3204.8072</pre>
=={{header|Julia}}==
<syntaxhighlight lang="julia"># Benford's Law
P(d) = log10(1+1/d)
function benford(numbers)
firstdigit(n) = last(digits(n))
counts = zeros(9)
foreach(n -> counts[firstdigit(n)] += 1, numbers)
counts ./ sum(counts)
end
struct Fibonacci end
Base.iterate(::Fibonacci, (a, b) = big.((0, 1))) = b, (b, a + b)
sample = Iterators.take(Fibonacci(), 1000)
observed = benford(sample) .* 100
expected = P.(1:9) .* 100
table = Real[1:9 observed expected]
using Plots
plot([observed expected]; title = "Benford's Law",
seriestype = [:bar :line], linewidth = [0 5],
xticks = 1:9, xlabel = "first digit", ylabel = "frequency %",
label = ["1000 Fibonacci numbers" "P(d)=log10(1+1/d)"])
using Printf
println("Benford's Law\nFrequency of first digit\nin 1000 Fibonacci numbers")
println("digit observed expected")
foreach(i -> @printf("%3d%9.2f%%%9.2f%%\n", table[i,:]...), 1:9)</syntaxhighlight>
{{Out}}
<pre>Benford's Law
Frequency of first digit
in 1000 Fibonacci numbers
digit observed expected
1 30.10% 30.10%
8 5.30% 5.12%
9 4.50% 4.58%</pre>
=={{header|Kotlin}}==
<syntaxhighlight lang="scala">import java.math.BigInteger
interface NumberGenerator {
Line 1,300 ⟶ 2,416:
}
class Benford(
override fun toString() = str
private val firstDigits = IntArray(9)
private val count = ng.numbers.size.toDouble()
private val str: String
init {
for (n in ng.numbers) {
firstDigits[n.toString().substring(0, 1).toInt() - 1]++
}
str = with(StringBuilder()) {
for (i in firstDigits.indices) {
append(i + 1).append('\t').append(firstDigits[i] / count)
append('\t').append(Math.log10(1 + 1.0 / (i + 1))).append('\n')
}
toString()
}
}
}
Line 1,322 ⟶ 2,442:
override val numbers: Array<BigInteger> by lazy {
val fib = Array<BigInteger>(1000, { BigInteger.ONE })
for (i in 2
fib[i] = fib[i - 2].add(fib[i - 1])
fib
Line 1,328 ⟶ 2,448:
}
fun main(a: Array<String>) = println(Benford(FibonacciGenerator))</
=={{header|Liberty BASIC}}==
Using function from
http://rosettacode.org/wiki/Fibonacci_sequence#Liberty_BASIC
<syntaxhighlight lang="lb">
dim bin(9)
Line 1,369 ⟶ 2,488:
fiboI = a
end function
</syntaxhighlight>
{{out}}
Line 1,384 ⟶ 2,503:
9 0.045 0.046
</pre>
=={{header|Lua}}==
<
expected = {}
for i = 1, 9 do
Line 1,405 ⟶ 2,523:
for i = 1, 9 do
print(i, actual[i] / n, expected[i])
end</
{{out}}
<pre>digit actual expected
Line 1,417 ⟶ 2,535:
8 0.053 0.051152522447381
9 0.045 0.045757490560675</pre>
=={{header|Mathematica}} / {{header|Wolfram Language}}==
<
Table[{d, N@Count[fibdata, d]/Length@fibdata, Log10[1. + 1/d]}, {d, 1,
9}] // Grid</
{{out}}
<pre>1 0.301 0.30103
Line 1,432 ⟶ 2,549:
8 0.053 0.0511525
9 0.045 0.0457575</pre>
=={{header|MATLAB}}==
{{trans|Julia}}
<syntaxhighlight lang="MATLAB">
benfords_law();
function benfords_law
% Benford's Law
P = @(d) log10(1 + 1./d);
% Benford function
function counts = benford(numbers)
firstdigit = @(n) floor(mod(n / 10^floor(log10(n)), 10));
counts = zeros(1, 9);
for i = 1:length(numbers)
digit = firstdigit(numbers(i));
if digit ~= 0
counts(digit) = counts(digit) + 1;
end
end
counts = counts ./ sum(counts);
end
% Generate Fibonacci numbers
function fibNums = fibonacci(n)
fibNums = zeros(1, n);
a = 0;
b = 1;
for i = 1:n
c = b;
b = a + b;
a = c;
fibNums(i) = b;
end
end
% Sample
sample = fibonacci(1000);
% Observed and expected frequencies
observed = benford(sample) * 100;
expected = arrayfun(P, 1:9) * 100;
% Table
mytable = [1:9; observed; expected]';
% Plotting
bar(1:9, observed);
hold on;
plot(1:9, expected, 'LineWidth', 2);
hold off;
title("Benford's Law");
xlabel("First Digit");
ylabel("Frequency %");
legend("1000 Fibonacci Numbers", "P(d) = log10(1 + 1/d)");
xticks(1:9);
% Displaying the results
fprintf("Benford's Law\nFrequency of first digit\nin 1000 Fibonacci numbers\n");
disp(table(mytable(:,1),mytable(:,2),mytable(:,3),'VariableNames',{'digit', 'observed(%)', 'expected(%)'}))
end
</syntaxhighlight>
{{out}}
<pre>
Benford's Law
Frequency of first digit
in 1000 Fibonacci numbers
digit observed(%) expected(%)
_____ ___________ ___________
1 30 30.103
2 17.7 17.609
3 12.5 12.494
4 9.6 9.691
5 8 7.9181
6 6.7 6.6947
7 5.7 5.7992
8 5.3 5.1153
9 4.5 4.5757
</pre>
=={{header|NetRexx}}==
<
options replace format comments java crossref symbols nobinary
Line 1,475 ⟶ 2,673:
brenfordDeveation(fibList)
return
</syntaxhighlight>
{{out}}
<pre>
Line 1,490 ⟶ 2,688:
9: 4.500000% 4.575749% 0.075749%
</pre>
=={{header|Nim}}==
<syntaxhighlight lang="nim">import math
import strformat
type
# Non zero digit range.
Digit = range[1..9]
# Count array used to compute a distribution.
Count = array[Digit, Natural]
# Array to store frequencies.
Distribution = array[Digit, float]
####################################################################################################
# Fibonacci numbers generation.
import bignum
proc fib(count: int): seq[Int] =
## Build a sequence of "count auccessive Finonacci numbers.
result.setLen(count)
result[0..1] = @[newInt(1), newInt(1)]
for i in 2..<count:
result[i] = result[i-1] + result[i-2]
####################################################################################################
# Benford distribution.
proc benfordDistrib(): Distribution =
## Compute Benford distribution.
for d in 1..9:
result[d] = log10(1 + 1 / d)
const BenfordDist = benfordDistrib()
#---------------------------------------------------------------------------------------------------
template firstDigit(n: Int): Digit =
# Return the first (non null) digit of "n".
ord(($n)[0]) - ord('0')
#---------------------------------------------------------------------------------------------------
proc actualDistrib(s: openArray[Int]): Distribution =
## Compute actual distribution of first digit.
## Null values are ignored.
var counts: Count
for val in s:
if not val.isZero():
inc counts[val.firstDigit()]
let total = sum(counts)
for d in 1..9:
result[d] = counts[d] / total
#---------------------------------------------------------------------------------------------------
proc display(distrib: Distribution) =
## Display an actual distribution versus the Benford reference distribution.
echo "d actual expected"
echo "---------------------"
for d in 1..9:
echo fmt"{d} {distrib[d]:6.4f} {BenfordDist[d]:6.4f}"
#———————————————————————————————————————————————————————————————————————————————————————————————————
let fibSeq = fib(1000)
let distrib = actualDistrib(fibSeq)
echo "Fibonacci numbers first digit distribution:\n"
distrib.display()</syntaxhighlight>
{{out}}
<pre>Fibonacci numbers first digit distribution:
d actual expected
---------------------
1 0.3010 0.3010
2 0.1770 0.1761
3 0.1250 0.1249
4 0.0960 0.0969
5 0.0800 0.0792
6 0.0670 0.0669
7 0.0560 0.0580
8 0.0530 0.0512
9 0.0450 0.0458</pre>
=={{header|Oberon-2}}==
{{Works with|oo2c version 2}}
<
MODULE BenfordLaw;
IMPORT
Line 1,537 ⟶ 2,827:
END
END BenfordLaw.
</syntaxhighlight>
{{out}}
<pre>
Line 1,552 ⟶ 2,842:
9 0.045 0.046
</pre>
=={{header|OCaml}}==
For the Fibonacci sequence, we use the function from
https://rosettacode.org/wiki/Fibonacci_sequence#Arbitrary_Precision<br>
Note the remark about the compilation of the program there.
<
open Num
Line 1,588 ⟶ 2,877:
List.iter (Printf.printf "%f ") (List.map benfords_law xvalues) ;
Printf.printf "\n" ;;
</syntaxhighlight>
{{out}}
<pre>
Line 1,595 ⟶ 2,884:
Prediction of Benford's law:
0.301030 0.176091 0.124939 0.096910 0.079181 0.066947 0.057992 0.051153 0.045757
</pre>
=={{header|PARI/GP}}==
<
my(t=vector(9,n,sum(i=1,#v,v[i]==n)));
print("Digit\tActual\tExpected");
Line 1,606 ⟶ 2,894:
lucas(n)=fibonacci(n-1)+fibonacci(n+1);
dist(fibonacci)
dist(lucas)</
{{out}}
<pre>Digit Actual Expected
Line 1,629 ⟶ 2,917:
8 51 51
9 46 46</pre>
=={{header|Pascal}}==
<syntaxhighlight lang="pascal">program fibFirstdigit;
{$IFDEF FPC}{$MODE Delphi}{$ELSE}{$APPTYPE CONSOLE}{$ENDIF}
uses
sysutils;
type
tDigitCount = array[0..9] of LongInt;
var
s: Ansistring;
dgtCnt,
expectedCnt : tDigitCount;
procedure GetFirstDigitFibonacci(var dgtCnt:tDigitCount;n:LongInt=1000);
//summing up only the first 9 digits
//n = 1000 -> difference to first 9 digits complete fib < 100 == 2 digits
var
a,b,c : LongWord;//about 9.6 decimals
Begin
for a in dgtCnt do dgtCnt[a] := 0;
a := 0;b := 1;
while n > 0 do
Begin
c := a+b;
//overflow? round and divide by base 10
IF c < a then
Begin a := (a+5) div 10;b := (b+5) div 10;c := a+b;end;
a := b;b := c;
s := IntToStr(a);inc(dgtCnt[Ord(s[1])-Ord('0')]);
dec(n);
end;
end;
procedure InitExpected(var dgtCnt:tDigitCount;n:LongInt=1000);
var
i: integer;
begin
for i := 1 to 9 do
dgtCnt[i] := trunc(n*ln(1 + 1 / i)/ln(10));
end;
var
reldiff: double;
i,cnt: integer;
begin
cnt := 1000;
InitExpected(expectedCnt,cnt);
GetFirstDigitFibonacci(dgtCnt,cnt);
writeln('Digit count expected rel diff');
For i := 1 to 9 do
Begin
reldiff := 100*(expectedCnt[i]-dgtCnt[i])/expectedCnt[i];
writeln(i:5,dgtCnt[i]:7,expectedCnt[i]:10,reldiff:10:5,' %');
end;
end.</syntaxhighlight>
<pre>Digit Count Expected rel Diff
1 301 301 0.00000 %
2 177 176 -0.56818 %
3 125 124 -0.80645 %
4 96 96 0.00000 %
5 80 79 -1.26582 %
6 67 66 -1.51515 %
7 56 57 1.75439 %
8 53 51 -3.92157 %
9 45 45 0.00000 %</pre>
=={{header|Perl}}==
<
use strict ;
use warnings ;
Line 1,655 ⟶ 3,007:
$result = sprintf ( "%.2f" , 100 * $expected[ $i - 1 ] ) ;
printf "%15s %%\n" , $result ;
}</
{{Out}}
<pre>
Line 1,669 ⟶ 3,021:
9 : 4.50 % 4.58 %
</pre>
=={{header|Phix}}==
{{trans|Go}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">benford</span><span style="color: #0000FF;">(</span><span style="color: #004080;">sequence</span> <span style="color: #000000;">s</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">string</span> <span style="color: #000000;">title</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">f</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">9</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">fdx</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sprint</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">])[</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]-</span><span style="color: #008000;">'0'</span>
<span style="color: #000000;">f</span><span style="color: #0000FF;">[</span><span style="color: #000000;">fdx</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">title</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Digit Observed% Predicted%"</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">o</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">]/</span><span style="color: #7060A8;">length</span><span style="color: #0000FF;">(</span><span style="color: #000000;">s</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">100</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">e</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">log10</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">/</span><span style="color: #000000;">i</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">100</span>
<span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">append</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">" %d %9.3f %8.3f"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">i</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">o</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">e</span><span style="color: #0000FF;">}))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">fib</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">lim</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">lim</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">lim</span> <span style="color: #008080;">do</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">res</span><span style="color: #0000FF;">[</span><span style="color: #000000;">i</span><span style="color: #0000FF;">],</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">+</span><span style="color: #000000;">a</span><span style="color: #0000FF;">}</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">res</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">res</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">benford</span><span style="color: #0000FF;">(</span><span style="color: #000000;">fib</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1000</span><span style="color: #0000FF;">),</span><span style="color: #008000;">"First 1000 Fibonacci numbers"</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">benford</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">get_primes</span><span style="color: #0000FF;">(-</span><span style="color: #000000;">10000</span><span style="color: #0000FF;">),</span><span style="color: #008000;">"First 10000 Prime numbers"</span><span style="color: #0000FF;">),</span>
<span style="color: #000000;">benford</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">sq_power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">3</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">tagset</span><span style="color: #0000FF;">(</span><span style="color: #000000;">500</span><span style="color: #0000FF;">)),</span><span style="color: #008000;">"First 500 powers of three"</span><span style="color: #0000FF;">)}</span>
<span style="color: #7060A8;">papply</span><span style="color: #0000FF;">(</span><span style="color: #004600;">true</span><span style="color: #0000FF;">,</span><span style="color: #7060A8;">printf</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,{</span><span style="color: #008000;">"%-40s%-40s%-40s\n"</span><span style="color: #0000FF;">},</span><span style="color: #7060A8;">columnize</span><span style="color: #0000FF;">(</span><span style="color: #000000;">res</span><span style="color: #0000FF;">)})</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
First 1000 Fibonacci numbers First 10000 Prime numbers First 500 powers of three
Line 1,781 ⟶ 3,068:
9 4.500 4.576 9 10.060 4.576 9 4.600 4.576
</pre>
=={{header|Picat}}==
<syntaxhighlight lang="picat">go =>
N = 1000,
Fib = [fib(I) : I in 1..N],
check_benford(Fib),
nl.
% Check a list of numbers for Benford's law
check_benford(L) =>
Len = L.len,
println(len=Len),
Count = [F[1].to_integer() : Num in L, F=Num.to_string()].occurrences(),
P = new_map([I=D/Len : I=D in Count]),
println("Benford (percent):"),
foreach(D in 1..9)
B = benford(D)*100,
PI = P.get(D,0)*100,
Diff = abs(PI - B),
printf("%d: count=%5d observed: %0.2f%% Benford: %0.2f%% diff=%0.3f\n", D,Count.get(D,0),PI,B,Diff)
end,
nl.
benford(D) = log10(1+1/D).
% create an occurrences map of a list
occurrences(List) = Map =>
Map = new_map(),
foreach(E in List)
Map.put(E, cond(Map.has_key(E),Map.get(E)+1,1))
end.</syntaxhighlight>
{{out}}
<pre>Benford (percent):
1: count= 301 observed: 30.10% Benford: 30.10% diff=0.003
2: count= 177 observed: 17.70% Benford: 17.61% diff=0.091
3: count= 125 observed: 12.50% Benford: 12.49% diff=0.006
4: count= 96 observed: 9.60% Benford: 9.69% diff=0.091
5: count= 80 observed: 8.00% Benford: 7.92% diff=0.082
6: count= 67 observed: 6.70% Benford: 6.69% diff=0.005
7: count= 56 observed: 5.60% Benford: 5.80% diff=0.199
8: count= 53 observed: 5.30% Benford: 5.12% diff=0.185
9: count= 45 observed: 4.50% Benford: 4.58% diff=0.076
</pre>
'''Extra credit:'''
Using data from https://en.wikipedia.org/wiki/Land_use_statistics_by_country "Cultivated land (km^2)"
{{out}}
<pre>Benford (percent):
1: count= 72 observed: 33.64% Benford: 30.10% diff=3.542
2: count= 38 observed: 17.76% Benford: 17.61% diff=0.148
3: count= 29 observed: 13.55% Benford: 12.49% diff=1.058
4: count= 19 observed: 8.88% Benford: 9.69% diff=0.812
5: count= 20 observed: 9.35% Benford: 7.92% diff=1.428
6: count= 12 observed: 5.61% Benford: 6.69% diff=1.087
7: count= 9 observed: 4.21% Benford: 5.80% diff=1.594
8: count= 11 observed: 5.14% Benford: 5.12% diff=0.025
9: count= 4 observed: 1.87% Benford: 4.58% diff=2.707</pre>
=={{header|PicoLisp}}==
Picolisp does not support floating point math, but it can call libc math functions and convert the results to a fixed point number for e.g. natural logarithm.
<syntaxhighlight lang="picolisp">
(scl 4)
(load "@lib/misc.l")
(load "@lib/math.l")
(setq LOG10E 0.4343)
(de fibo (N)
(let (A 0 B 1 C NIL)
(make
(link B)
(do (dec N)
(setq C (+ A B) A B B C)
(link B)))))
(setq Actual
(let (
Digits (sort (mapcar '((N) (format (car (chop N)))) (fibo 1000)))
Count 0
)
(make
(for (Ds Digits Ds (cdr Ds))
(inc 'Count)
(when (n== (car Ds) (cadr Ds))
(link Count)
(setq Count 0))))))
(setq Expected
(mapcar
'((D) (*/ (log (+ 1. (/ 1. D))) LOG10E 1.))
(range 1 9)))
(prinl "Digit\tActual\tExpected")
(let (As Actual Bs Expected)
(for D 9
(prinl D "\t" (format (pop 'As) 3) "\t" (round (pop 'Bs) 3))))
(bye)
</syntaxhighlight>
{{Out}}
<pre>
Digit Actual Expected
1 0.301 0.301
2 0.177 0.176
3 0.125 0.125
4 0.096 0.097
5 0.080 0.079
6 0.067 0.067
7 0.056 0.058
8 0.053 0.051
9 0.045 0.046
</pre>
=={{header|PL/I}}==
<syntaxhighlight lang="pl/i">
(fofl, size, subrg):
Benford: procedure options(main); /* 20 October 2013 */
Line 1,823 ⟶ 3,223:
end tally;
end Benford;
</syntaxhighlight>
Results:
<pre>
Line 1,837 ⟶ 3,237:
9 0.04575749 0.04499817
</pre>
=={{header|PL/pgSQL}}==
<syntaxhighlight lang="sql">
WITH recursive
constant(val) AS
Line 1,874 ⟶ 3,273:
(select cast(corr(probability_theoretical,probability_real) as numeric(5,4)) correlation
from benford) c
</syntaxhighlight>
=={{header|PowerShell}}==
The sample file was not found. I selected another that contained the first two-thousand in the Fibonacci sequence, so there is a small amount of extra filtering.
<syntaxhighlight lang="powershell">
$url = "https://oeis.org/A000045/b000045.txt"
$file = "$env:TEMP\FibonacciNumbers.txt"
Line 1,895 ⟶ 3,293:
Remove-Item -Path $file -Force -ErrorAction SilentlyContinue
</syntaxhighlight>
{{Out}}
<pre>
Line 1,910 ⟶ 3,308:
9 45 0.045 0.04576
</pre>
=={{header|Prolog}}==
{{works with|SWI Prolog|6.2.6 by Jan Wielemaker, University of Amsterdam}}
Note: SWI Prolog implements arbitrary precision integer arithmetic through use of the GNU MP library
<
% Does the Fibonacci sequence follow Benford's law?
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Line 1,957 ⟶ 3,354:
findall(B, (between(1,9,N), benford(N,B)), Benford),
findall(C, firstchar(C), Fc), freq(Fc, Freq),
writeHdr, maplist(writeData, Benford, Freq).</
{{out}}
<pre>?- go.
Line 1,970 ⟶ 3,367:
5.12% - 5.30%
4.58% - 4.50%</pre>
=={{header|PureBasic}}==
<syntaxhighlight lang="purebasic">#MAX_N=1000
NewMap d1.i()
Dim fi.s(#MAX_N)
fi(0)="0" : fi(1)="1"
Declare.s Sigma(sx.s,sy.s)
For I=2 To #MAX_N
fi(I)=Sigma(fi(I-2),fi(I-1))
Next
For I=1 To #MAX_N
d1(Left(fi(I),1))+1
Next
Procedure.s Sigma(sx.s, sy.s)
Define i.i, v1.i, v2.i, r.i
Define s.s, sa.s
sy=ReverseString(sy) : s=ReverseString(sx)
For i=1 To Len(s)*Bool(Len(s)>Len(sy))+Len(sy)*Bool(Len(sy)>=Len(s))
v1=Val(Mid(s,i,1))
v2=Val(Mid(sy,i,1))
r+v1+v2
sa+Str(r%10)
r/10
Next i
If r : sa+Str(r%10) : EndIf
ProcedureReturn ReverseString(sa)
EndProcedure
OpenConsole("Benford's law: Fibonacci sequence 1.."+Str(#MAX_N))
Print(~"Dig.\t\tCnt."+~"\t\tExp.\t\tDif.\n\n")
ForEach d1()
Print(RSet(MapKey(d1()),4," ")+~"\t:\t"+RSet(Str(d1()),3," ")+~"\t\t")
ex=Int(#MAX_N*Log(1+1/Val(MapKey(d1())))/Log(10))
PrintN(RSet(Str(ex),3," ")+~"\t\t"+RSet(StrF((ex-d1())*100/ex,5),8," ")+" %")
Next
PrintN(~"\nPress Enter...")
Input()</syntaxhighlight>
{{out}}
<pre>Dig. Cnt. Exp. Dif.
1 : 301 301 0.00000 %
2 : 177 176 -0.56818 %
3 : 125 124 -0.80645 %
4 : 96 96 0.00000 %
5 : 80 79 -1.26582 %
6 : 67 66 -1.51515 %
7 : 56 57 1.75439 %
8 : 53 51 -3.92157 %
9 : 45 45 0.00000 %
Press Enter...</pre>
=={{header|Python}}==
Works with Python 3.X & 2.7
<
from itertools import islice, count
from collections import Counter
Line 2,011 ⟶ 3,462:
# just to show that not all kind-of-random sets behave like that
show_dist("random", islice(heads(rand1000()), 10000))</
{{out}}
<pre>fibbed Benfords deviation
Line 2,045 ⟶ 3,496:
11.0% 5.1% 5.9%
10.9% 4.6% 6.3%</pre>
=={{header|R}}==
<syntaxhighlight lang="r">
pbenford <- function(d){
return(log10(1+(1/d)))
Line 2,079 ⟶ 3,529:
print(data)
</syntaxhighlight>
{{out}}
digit obs.frequency exp.frequency dev_percentage
Line 2,091 ⟶ 3,541:
8 0.053 0.05115 0.184748
9 0.045 0.04576 0.075749
=={{header|Racket}}==
<
(define (log10 n) (/ (log n) (log 10)))
Line 2,128 ⟶ 3,577:
;; 7: 5.8% 5.8%
;; 8: 5.1% 5.1%
;; 9: 4.6% 4.6%</
=={{header|Raku}}==
(formerly Perl 6)
{{Works with|rakudo|2016-10-24}}
<syntaxhighlight lang="raku" line>sub benford(@a) { bag +« @a».substr(0,1) }
sub show(%distribution) {
printf "%9s %9s %s\n", <Actual Expected Deviation>;
for 1 .. 9 -> $digit {
my $actual = %distribution{$digit} * 100 / [+] %distribution.values;
my $expected = (1 + 1 / $digit).log(10) * 100;
printf "%d: %5.2f%% | %5.2f%% | %.2f%%\n",
$digit, $actual, $expected, abs($expected - $actual);
}
}
multi MAIN($file) { show benford $file.IO.lines }
multi MAIN() { show benford ( 1, 1, 2, *+* ... * )[^1000] }</syntaxhighlight>
'''Output:''' First 1000 Fibonaccis
<pre> Actual Expected Deviation
1: 30.10% | 30.10% | 0.00%
2: 17.70% | 17.61% | 0.09%
3: 12.50% | 12.49% | 0.01%
4: 9.60% | 9.69% | 0.09%
5: 8.00% | 7.92% | 0.08%
6: 6.70% | 6.69% | 0.01%
7: 5.60% | 5.80% | 0.20%
8: 5.30% | 5.12% | 0.18%
9: 4.50% | 4.58% | 0.08%</pre>
'''Extra credit:''' Square Kilometers of land under cultivation, by country / territory. First column from Wikipedia: [[wp:Land_use_statistics_by_country|Land use statistics by country]].
<pre> Actual Expected Deviation
1: 33.33% | 30.10% | 3.23%
2: 18.31% | 17.61% | 0.70%
3: 13.15% | 12.49% | 0.65%
4: 8.45% | 9.69% | 1.24%
5: 9.39% | 7.92% | 1.47%
6: 5.63% | 6.69% | 1.06%
7: 4.69% | 5.80% | 1.10%
8: 5.16% | 5.12% | 0.05%
9: 1.88% | 4.58% | 2.70%</pre>
=={{header|REXX}}==
The REXX language (for the most part) hasn't any high math functions, so the '''e''', '''ln''', and '''log''' functions were included herein.
For the extra credit stuff, it was chosen to generate Fibonacci and factorials rather than find a
web─page with them listed, as each list is very easy to generate.
<syntaxhighlight lang="rexx">/*REXX pgm demonstrates Benford's law applied to 2 common functions (30 dec. digs used).*/
numeric digits length( e() ) - length(.) /*width of (e) for LN & LOG precision.*/
parse arg N .; if N=='' | N=="," then N= 1000 /*allow sample size to be specified. */
pad= " " /*W1, W2: # digs past the decimal point*/
w1= max(2 + length('observed'), length(N-2) ) /*for aligning output for a number. */
w2= max(2 + length('expected'), length(N ) ) /* " " frequency distributions.*/
LN10= ln(10) /*calculate the ln(10) {used by LOG}*/
call coef /*generate nine frequency coefficients.*/
call fib /*generate N Fibonacci numbers. */
call show "Benford's law applied to" N 'Fibonacci numbers'
call fact /*generate N factorials. */
call show "Benford's law applied to" N 'factorial products'
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
fact: @.=1; do j=2 for N-1; a= j-1; @.j= @.a * j; end; return
fib:
e: return 2.71828182845904523536028747135266249775724709369995957496696762772407663035
log: return ln( arg(1) ) / LN10
/*──────────────────────────────────────────────────────────────────────────────────────*/
ln: procedure; parse arg x; e= e(); _= e; ig= (x>1.5); is= 1 - 2 * (ig\=1); i= 0; s= x
do while ig&s>1.5 | \ig&s<.5 /*nitty─gritty part of LN calculation*/
do k=-1; iz=s*_**-is; if k>=0&(ig&iz<1|\ig&iz>.5) then leave; _=_*_; izz=iz; end
s=izz; i= i + is* 2**k; end /*while*/; x= x * e** - i - 1; z= 0; _= -1; p= z
do k=1; _= -_ * x; z= z + _/k; if z=p then leave; p= z; end /*k*/; return z+i
/*──────────────────────────────────────────────────────────────────────────────────────*/
show: say; say pad ' digit ' pad center("observed",w1) pad center('expected',w2)
say pad '───────' pad center("", w1, '─') pad center("",w2,'─') pad arg(1)
!.=0; do j=1 for N; _= left(@.j, 1); !._= !._ + 1 /*get the 1st digit.*/
end /*j*/
do f=1 for 9; say pad center(f,7) pad center(format(!.f/N,,length(N-2)),w1) #.f
end /*k*/
return</syntaxhighlight>
{{out|output|text= when using the default (1000 numbers) for the input:}}
<pre>
digit observed expected
─────── ────────── ────────── Benford's law applied to 1000 Fibonacci numbers
1 0.301 0.301030
2 0.177 0.176091
3 0.125 0.124939
4 0.096 0.096910
5 0.080 0.079181
6 0.067 0.066947
7 0.056 0.057992
8 0.053 0.051153
9 0.045 0.045757
digit observed expected
─────── ────────── ────────── Benford's law applied to 1000 factorial products
1 0.293
2 0.176
3
4 0.102 0.096910
5 0.069 0.079181
6 0.087 0.066947
7 0.051 0.057992
8 0.051 0.051153
9 0.047 0.045757
</pre>
{{out|output|text= when using the following for the input: <tt> 10000 </tt>}}
<pre>
digit observed expected
1 0.
2 0.
3 0.
4 0.
5 0.
6 0.
7 0.
8 0.
9 0.
digit observed expected
1 0.
2 0.
3 0.
4 0.
5 0.
6 0.
7 0.
8 0.
9 0.
</pre>
=={{header|Ring}}==
<syntaxhighlight lang="ring">
# Project : Benford's law
decimals(3)
n= 1000
actual = list(n)
for x = 1 to len(actual)
actual[x] = 0
next
for nr = 1 to n
n1 = string(fibonacci(nr))
j = number(left(n1,1))
actual[j] = actual[j] + 1
next
see "Digit " + "Actual " + "Expected" + nl
for m = 1 to 9
fr = frequency(m)*100
see "" + m + " " + (actual[m]/10) + " " + fr + nl
next
func frequency(n)
freq = log10(n+1) - log10(n)
return freq
func log10(n)
log1 = log(n) / log(10)
return log1
func fibonacci(y)
if y = 0 return 0 ok
if y = 1 return 1 ok
if y > 1 return fibonacci(y-1) + fibonacci(y-2) ok
</syntaxhighlight>
Output:
<pre>
Digit Actual Expected
1 30.100 30.103
2 17.700 17.609
3 12.500 12.494
4 9.500 9.691
5 8.000 7.918
6 6.700 6.695
7 5.600 5.799
8 5.300 5.115
9 4.500 4.576
</pre>
=={{header|RPL}}==
Here we use an interesting RPL feature that allows to pass an arithmetic expression or a function name as an argument. Thus, the code below can check the 'benfordance' of various series without needing to edit it. The resulting array being a matrix, some additional operations allow to compute the maximum difference at digit level between the tested function and Benford's law, as a distance indicator.
{{works with|Halcyon Calc|4.2.7}}
≪ → sn min max
≪ { 9 2 } 0 CON
min max FOR j
{ 0 1 } 1 j sn EVAL
IF DUP THEN
MANT FLOOR PUT
DUP2 GET 1 + PUT
ELSE
3 DROPN
END
NEXT
max min - 1 + /
1 9 FOR j
{ 0 2 } 1 j PUT
j INV 1 + LOG PUT
NEXT
DUP [ 1 0 ] *
OVER [ 0 1 ] * -
RNRM
≫
≫
'BENFD' STO
'→FIB' 1 100 BENFD
≪ LOG ≫ 1 10000 BENFD
Output:
<pre>
4: [[ 0.301 0.301029995664 ]
[ 0.177 0.176091259056 ]
[
[ 6.7E-02 6.69467896306E-02 ]
[
[ 0.045 4.57574905607E-02 ]]
3: 1.8847477552E-02
2: [[ 9E-04 0.301029995664 ]
[ 9E-03 0.176091259056 ]
[ 0.09001 0.124938736608 ]
[ 0.90001 9.69100130081E-02 ]
[ 0.00001 7.91812460476E-02 ]
[ 0.00002 6.69467896306E-02 ]
[ 0.00001 5.79919469777E-02 ]
[ 0.00001 5.11525224474E-02 ]
[ 0.00002 4.57574905607E-02 ]]
1: 0.775161263392
</pre>
=={{header|Ruby}}==
{{trans|Python}}
<
def fib(n)
Line 2,265 ⟶ 3,850:
# just to show that not all kind-of-random sets behave like that
show_dist("random", random(10000))</
{{out}}
Line 2,302 ⟶ 3,887:
9: 10.8% 4.6% 6.2%
</pre>
=={{header|Run BASIC}}==
<
N = 1000
for i = 0 to N - 1
Line 2,335 ⟶ 3,919:
next i
end function
</syntaxhighlight>
<table border=1><TR bgcolor=wheat><TD>Digit<td>Actual<td>Expected</td><tr><tr align=right><td>1</td><td>30.100</td><td>30.103</td></tr><tr align=right><td>2</td><td>17.700</td><td>17.609</td></tr><tr align=right><td>3</td><td>12.500</td><td>12.494</td></tr><tr align=right><td>4</td><td> 9.500</td><td> 9.691</td></tr><tr align=right><td>5</td><td> 8.000</td><td> 7.918</td></tr><tr align=right><td>6</td><td> 6.700</td><td> 6.695</td></tr><tr align=right><td>7</td><td> 5.600</td><td> 5.799</td></tr><tr align=right><td>8</td><td> 5.300</td><td> 5.115</td></tr><tr align=right><td>9</td><td> 4.500</td><td> 4.576</td></tr></table>
=={{header|Rust}}==
{{works with|rustc|1.12 stable}}
This solution uses the ''num'' create for arbitrary-precision integers and the ''num_traits'' create for the ''zero'' and ''one'' implementations. It computes the Fibonacci numbers from scratch via the ''fib'' function.
<syntaxhighlight lang="rust">
extern crate num_traits;
extern crate num;
use num::bigint::{BigInt, ToBigInt};
use num_traits::{Zero, One};
use std::collections::HashMap;
// Return a vector of all fibonacci results from fib(1) to fib(n)
fn fib(n: usize) -> Vec<BigInt> {
let mut result = Vec::with_capacity(n);
let mut a = BigInt::zero();
let mut b = BigInt::one();
result.push(b.clone());
for i in 1..n {
let t = b.clone();
b = a+b;
a = t;
result.push(b.clone());
}
result
}
// Return the first digit of a `BigInt`
fn first_digit(x: &BigInt) -> u8 {
let zero = BigInt::zero();
assert!(x > &zero);
let s = x.to_str_radix(10);
// parse the first digit of the stringified integer
*&s[..1].parse::<u8>().unwrap()
}
fn main() {
const N: usize = 1000;
let mut counter: HashMap<u8, u32> = HashMap::new();
for x in fib(N) {
let d = first_digit(&x);
*counter.entry(d).or_insert(0) += 1;
}
println!("{:>13} {:>10}", "real", "predicted");
for y in 1..10 {
println!("{}: {:10.3} v. {:10.3}", y, *counter.get(&y).unwrap_or(&0) as f32 / N as f32,
(1.0 + 1.0 / (y as f32)).log10());
}
}
</syntaxhighlight>
{{out}}
<pre>
real predicted
1: 0.301 v. 0.301
2: 0.177 v. 0.176
3: 0.125 v. 0.125
4: 0.096 v. 0.097
5: 0.080 v. 0.079
6: 0.067 v. 0.067
7: 0.056 v. 0.058
8: 0.053 v. 0.051
9: 0.045 v. 0.046
</pre>
=={{header|Scala}}==
<
val fibs : Stream[BigInt] = { def series(i:BigInt,j:BigInt):Stream[BigInt] = i #:: series(j, i+j); series(1,0).tail.tail }
Line 2,367 ⟶ 4,021:
case (k, v) => println( "%d: %5.2f%% | %5.2f%% | %5.4f%%".format(k,v._2*100,v._1*100,math.abs(v._2-v._1)*100) )
}
}</
{{out}}
<pre>
Line 2,383 ⟶ 4,037:
9: 4.50% | 4.58% | 0.0757%
</pre>
=={{header|Scheme}}==
{{works with|Chez Scheme}}
<syntaxhighlight lang="scheme">; Compute the probability of leading digit d (an integer [1,9]) according to Benford's law.
(define benford-probability
(lambda (d)
(- (log (1+ d) 10) (log d 10))))
; Determine the distribution of the leading digit of the sequence provided by the given
; generator. Return as a vector of 10 elements -- the 0th element will always be 0.
(define leading-digit-distribution
(lambda (seqgen count)
(let ((digcounts (make-vector 10 0)))
(do ((index 0 (1+ index)))
((>= index count))
(let* ((value (seqgen))
(string (format "~a" value))
(leadchr (string-ref string 0))
(leaddig (- (char->integer leadchr) (char->integer #\0))))
(vector-set! digcounts leaddig (1+ (vector-ref digcounts leaddig)))))
(vector-map (lambda (digcnt) (/ digcnt count)) digcounts))))
; Create a Fibonacci sequence generator.
(define make-fibgen
(lambda ()
(let ((fn-2 0) (fn-1 1))
(lambda ()
(let ((fn fn-1))
(set! fn-1 (+ fn-2 fn-1))
(set! fn-2 fn)
fn)))))
; Create a sequence generator that returns elements of the given vector.
(define make-vecgen
(lambda (vector)
(let ((index 0))
(lambda ()
(set! index (1+ index))
(vector-ref vector (1- index))))))
; Read all the values in the given file into a list.
(define list-read-file
(lambda (filenm)
(call-with-input-file filenm
(lambda (ip)
(let accrue ((value (read ip)))
(if (eof-object? value)
'()
(cons value (accrue (read ip)))))))))
; Display a table of digit, Benford's law, sequence distribution, and difference.
(define display-table
(lambda (seqnam seqgen count)
(printf "~%~3@a ~11@a ~11@a ~11@a~%" "dig" "Benford's" seqnam "difference")
(let ((dist (leading-digit-distribution seqgen count)))
(do ((digit 1 (1+ digit)))
((> digit 9))
(let* ((fraction (vector-ref dist digit))
(benford (benford-probability digit))
(diff (- fraction benford)))
(printf "~3d ~11,5f ~11,5f ~11,5f~%" digit benford fraction diff))))))
; Emit tables of various sequence distributions.
(display-table "Fib/1000" (make-fibgen) 1000)
(display-table "Rnd/1T/1M" (lambda () (1+ (random 1000000000000))) 1000000)
(let ((craters (list->vector (list-read-file "moon_craters.lst"))))
(display-table "Craters/D" (make-vecgen craters) (vector-length craters)))</syntaxhighlight>
{{out}}
The first one thousand Fibonnaci numbers.
<pre>dig Benford's Fib/1000 difference
1 0.30103 0.30100 -0.00003
2 0.17609 0.17700 0.00091
3 0.12494 0.12500 0.00006
4 0.09691 0.09600 -0.00091
5 0.07918 0.08000 0.00082
6 0.06695 0.06700 0.00005
7 0.05799 0.05600 -0.00199
8 0.05115 0.05300 0.00185
9 0.04576 0.04500 -0.00076</pre>
{{out}}
One million pseudo-random numbers from one to one trillion.
<pre>dig Benford's Rnd/1T/1M difference
1 0.30103 0.11055 -0.19048
2 0.17609 0.11100 -0.06509
3 0.12494 0.11126 -0.01368
4 0.09691 0.11129 0.01438
5 0.07918 0.11141 0.03223
6 0.06695 0.11088 0.04394
7 0.05799 0.11120 0.05321
8 0.05115 0.11112 0.05997
9 0.04576 0.11128 0.06553</pre>
{{out}}
Diameters of 1,595 Lunar craters in meters.
<br />
From the Gazetteer of Planetary Nomenclature [https://planetarynames.wr.usgs.gov/],
referred to by the Wikipedia Planetary nomenclature page [https://en.wikipedia.org/wiki/Planetary_nomenclature].
<pre>dig Benford's Craters/D difference
1 0.30103 0.23950 -0.06153
2 0.17609 0.10658 -0.06951
3 0.12494 0.11912 -0.00582
4 0.09691 0.12414 0.02723
5 0.07918 0.11034 0.03116
6 0.06695 0.10596 0.03901
7 0.05799 0.06959 0.01160
8 0.05115 0.06646 0.01531
9 0.04576 0.05831 0.01255</pre>
=={{header|SETL}}==
<syntaxhighlight lang="setl">program benfords_law;
fibos := fibo_list(1000);
expected := [log(1 + 1/d)/log 10 : d in [1..9]];
actual := benford(fibos);
print('d Expected Actual');
loop for d in [1..9] do
print(d, ' ', fixed(expected(d), 7, 5), ' ', fixed(actual(d), 7, 5));
end loop;
proc benford(list);
dist := [];
loop for n in list do
dist(val(str n)(1)) +:= 1;
end loop;
return [d / #list : d in dist];
end proc;
proc fibo_list(n);
a := 1;
b := 1;
fibs := [];
loop while n>0 do
fibs with:= a;
[a, b] := [b, a+b];
n -:= 1;
end loop;
return fibs;
end proc;
end program;</syntaxhighlight>
{{out}}
<pre>d Expected Actual
1 0.30103 0.30100
2 0.17609 0.17700
3 0.12494 0.12500
4 0.09691 0.09600
5 0.07918 0.08000
6 0.06695 0.06700
7 0.05799 0.05600
8 0.05115 0.05300
9 0.04576 0.04500</pre>
=={{header|Sidef}}==
<
var fibonacci = 1000.of {|i| fib(i).digit(-1) }
for i in (1..9) {
var num = fibonacci.count_by {|j| j == i }
actuals.append(num / 1000)
expected.append(1 + (1/i) -> log10)
}
"%17s%17s\n".printf("Observed","Expected")
for i in (1..9) {
"%d : %11s %%%15s %%\n".printf(
i, "%.2f".sprintf(100 * actuals[i - 1]),
"%.2f".sprintf(100 * expected[i - 1]),
)
}</syntaxhighlight>
{{out}}
<pre>
Line 2,426 ⟶ 4,230:
The query is the same for any number sequence you care to put in the <tt>benford</tt> table.
<
create table benford (num integer);
Line 2,468 ⟶ 4,272:
-- Tidy up
drop table benford;</
{{out}}
Line 2,485 ⟶ 4,289:
9 rows selected.</pre>
=={{header|Stata}}==
<syntaxhighlight lang="stata">clear
set obs 1000
scalar phi=(1+sqrt(5))/2
gen fib=(phi^_n-(-1/phi)^_n)/sqrt(5)
gen k=real(substr(string(fib),1,1))
hist k, discrete // show a histogram
qui tabulate k, matcell(f) // compute frequencies
mata
f=st_matrix("f")
p=log10(1:+1:/(1::9))*sum(f)
// print observed vs predicted probabilities
f,p
1 2
+-----------------------------+
1 | 297 301.0299957 |
2 | 178 176.0912591 |
3 | 127 124.9387366 |
4 | 96 96.91001301 |
5 | 80 79.18124605 |
6 | 67 66.94678963 |
7 | 57 57.99194698 |
8 | 53 51.15252245 |
9 | 45 45.75749056 |
+-----------------------------+</syntaxhighlight>
Assuming the data are random, one can also do a goodness of fit [https://en.wikipedia.org/wiki/Pearson%27s_chi-squared_test chi-square test]:
<syntaxhighlight lang="stata">// chi-square statistic
chisq=sum((f-p):^2:/p)
chisq
.2219340262
// p-value
chi2tail(8,chisq)
.9999942179
end</syntaxhighlight>
The p-value is very close to 1, showing that the observed distribution is very close to the Benford law.
The fit is not as good with the sequence (2+sqrt(2))^n:
<syntaxhighlight lang="stata">clear
set obs 500
scalar s=2+sqrt(2)
gen a=s^_n
gen k=real(substr(string(a),1,1))
hist k, discrete
qui tabulate k, matcell(f)
mata
f=st_matrix("f")
p=log10(1:+1:/(1::9))*sum(f)
f,p
1 2
+-----------------------------+
1 | 134 150.5149978 |
2 | 99 88.04562953 |
3 | 68 62.4693683 |
4 | 34 48.4550065 |
5 | 33 39.59062302 |
6 | 33 33.47339482 |
7 | 33 28.99597349 |
8 | 33 25.57626122 |
9 | 33 22.87874528 |
+-----------------------------+
chisq=sum((f-p):^2:/p)
chisq
16.26588528
chi2tail(8,chisq)
.0387287805
end</syntaxhighlight>
Now the p-value is less than the usual 5% risk, and one would reject the hypothesis that the data follow the Benford law.
=={{header|Swift}}==
<syntaxhighlight lang="swift">import Foundation
/* Reads from a file and returns the content as a String */
func readFromFile(fileName file:String) -> String{
var ret:String = ""
let path = Foundation.URL(string: "file://"+file)
do {
ret = try String(contentsOf: path!, encoding: String.Encoding.utf8)
}
catch {
print("Could not read from file!")
exit(-1)
}
return ret
}
/* Calculates the probability following Benford's law */
func benford(digit z:Int) -> Double {
if z<=0 || z>9 {
perror("Argument must be between 1 and 9.")
return 0
}
return log10(Double(1)+Double(1)/Double(z))
}
// get CLI input
if CommandLine.arguments.count < 2 {
print("Usage: Benford [FILE]")
exit(-1)
}
let pathToFile = CommandLine.arguments[1]
// Read from given file and parse into lines
let content = readFromFile(fileName: pathToFile)
let lines = content.components(separatedBy: "\n")
var digitCount:UInt64 = 0
var countDigit:[UInt64] = [0,0,0,0,0,0,0,0,0]
// check digits line by line
for line in lines {
if line == "" {
continue
}
let charLine = Array(line.characters)
switch(charLine[0]){
case "1":
countDigit[0] += 1
digitCount += 1
break
case "2":
countDigit[1] += 1
digitCount += 1
break
case "3":
countDigit[2] += 1
digitCount += 1
break
case "4":
countDigit[3] += 1
digitCount += 1
break
case "5":
countDigit[4] += 1
digitCount += 1
break
case "6":
countDigit[5] += 1
digitCount += 1
break
case "7":
countDigit[6] += 1
digitCount += 1
break
case "8":
countDigit[7] += 1
digitCount += 1
break
case "9":
countDigit[8] += 1
digitCount += 1
break
default:
break
}
}
// print result
print("Digit\tBenford [%]\tObserved [%]\tDeviation")
print("~~~~~\t~~~~~~~~~~~~\t~~~~~~~~~~~~\t~~~~~~~~~")
for i in 0..<9 {
let temp:Double = Double(countDigit[i])/Double(digitCount)
let ben = benford(digit: i+1)
print(String(format: "%d\t%.2f\t\t%.2f\t\t%.4f", i+1,ben*100,temp*100,ben-temp))
}</syntaxhighlight>
{{out}}
<pre>$ ./Benford
Usage: Benford [FILE]
$ ./Benford Fibonacci.txt
Digit Benford [%] Observed [%] Deviation
~~~~~ ~~~~~~~~~~~~ ~~~~~~~~~~~~ ~~~~~~~~~
1 30.10 30.10 0.0000
2 17.61 17.70 -0.0009
3 12.49 12.50 -0.0001
4 9.69 9.60 0.0009
5 7.92 8.00 -0.0008
6 6.69 6.70 -0.0001
7 5.80 5.60 0.0020
8 5.12 5.30 -0.0018
9 4.58 4.50 0.0008</pre>
=={{header|Tcl}}==
<
# Count the leading digits (RE matches first digit in each number,
# even if negative)
Line 2,505 ⟶ 4,505:
[expr {log(1+1./$digit)/log(10)*100.0}]]
}
}</
Demonstrating with Fibonacci numbers:
<
for {set a 1;set b [set i 0]} {$i < $n} {incr i} {
lappend result [set b [expr {$a + [set a $b]}]]
Line 2,513 ⟶ 4,513:
return $result
}
benfordTest [fibs 1000]</
{{out}}
<pre>
Line 2,529 ⟶ 4,529:
</pre>
=={{header|True BASIC}}==
{{trans|Chipmunk Basic}}
<syntaxhighlight lang="qbasic">FUNCTION log10(n)
LET log10 = LOG(n)/LOG(10)
END FUNCTION
FUNCTION frequency(n)
LET frequency = (log10(n+1)-log10(n))
END FUNCTION
FUNCTION fibonacci(n)
!https://rosettacode.org/wiki/Fibonacci_sequence#True_BASIC
LET n1 = 0
LET n2 = 1
FOR k = 1 TO ABS(n)
LET sum = n1+n2
LET n1 = n2
LET n2 = sum
NEXT k
IF n < 0 THEN LET fibonacci = n1*((-1)^((-n)+1)) ELSE LET fibonacci = n1
END FUNCTION
CLEAR
LET n = 1000
DIM actual(0)
MAT REDIM actual(n)
FOR nr = 1 TO n
LET num$ = STR$(fibonacci(nr))
LET j = VAL((num$)[1:1])
LET actual(j) = actual(j)+1
NEXT nr
PRINT "First 1000 Fibonacci numbers"
PRINT "Digit Actual freq Expected freq"
FOR i = 1 TO 9
LET freq = frequency(i)*100
PRINT USING "###": i;
PRINT USING " ##.###": actual(i)/10;
PRINT USING " ##.###": freq
NEXT i
END</syntaxhighlight>
=={{header|VBA (Visual Basic for Application)}}==
<syntaxhighlight lang="vb">
Sub BenfordLaw()
Dim BenResult(1 To 9) As Long
BENref = "30,1%|17,6%|12,5%|9,7%|7,9%|6,7%|5,8%|5,1%|4,6%"
For Each c In Selection.Cells
If InStr(1, "-0123456789", Left(c, 1)) > 0 Then
For i = 1 To 9
If CInt(Left(Abs(c), 1)) = i Then BenResult(i) = BenResult(i) + 1: Exit For
Next
End If
Next
Total= Application.Sum(BenResult)
biggest= Len(CStr(BenResult(1)))
txt = "# | Values | Real | Expected " & vbCrLf
For i = 1 To 9
If BenResult(i) > 0 Then
txt = txt & "#" & i & " | " & vbTab
txt = txt & String((biggest - Len(CStr(BenResult(i)))) * 2, " ") & Format(BenResult(i), "0") & " | " & vbTab
txt = txt & String((Len(CStr(Format(BenResult(1) / Total, "##0.0%"))) - Len(CStr(Format(BenResult(i) / Total, "##0.0%")))) * 2, " ") & Format(BenResult(i) / Total, "##0.0%") & " | " & vbTab
txt = txt & Format(Split(BENref, "|")(i - 1), " ##0.0%") & vbCrLf
End If
Next
MsgBox txt, vbOKOnly, "Finish"
End Sub
}
</syntaxhighlight>
=={{header|Visual FoxPro}}==
<
#DEFINE CTAB CHR(9)
#DEFINE COMMA ","
Line 2,571 ⟶ 4,646:
RETURN LOG10(1 + 1/d)
ENDFUNC
</syntaxhighlight>
{{out}}
<pre>
Line 2,587 ⟶ 4,662:
Correlation Coefficient: 0.999908
</pre>
=={{header|V (Vlang)}}==
{{trans|Go}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="v (vlang)">import math
fn fib1000() []f64 {
mut a, mut b, mut r := 0.0, 1.0, []f64{len:1000}
for i in 0..r.len {
r[i], a, b = b, b, b+a
}
return r
}
fn main() {
show(fib1000(), "First 1000 Fibonacci numbers")
}
fn show(c []f64, title string) {
mut f := [9]int{}
for v in c {
f["$v"[0]-'1'[0]]++
}
println(title)
println("Digit Observed Predicted")
for i, n in f {
println(" ${i+1} ${f64(n)/f64(c.len):9.3f} ${math.log10(1+1/f64(i+1)):8.3f}")
}
}</syntaxhighlight>
{{out}}
<pre>
First 1000 Fibonacci numbers:
Digit Observed Predicted
1 0.301 0.301
2 0.177 0.176
3 0.125 0.125
4 0.096 0.097
5 0.080 0.079
6 0.067 0.067
7 0.056 0.058
8 0.053 0.051
9 0.045 0.046
</pre>
=={{header|Wren}}==
{{trans|Go}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="wren">import "./fmt" for Fmt
var fib1000 = Fn.new {
var a = 0
var b = 1
var r = List.filled(1000, 0)
for (i in 0...r.count) {
var oa = a
var ob = b
r[i] = ob
a = ob
b = ob + oa
}
return r
}
var LN10 = 2.3025850929940457
var log10 = Fn.new { |x| x.log / LN10 }
var show = Fn.new { |c, title|
var f = List.filled(9, 0)
for (v in c) {
var t = "%(v)".bytes[0] - 49
f[t] = f[t] + 1
}
System.print(title)
System.print("Digit Observed Predicted")
for (i in 0...f.count) {
var n = f[i]
var obs = Fmt.f(9, n/c.count, 3)
var t = log10.call(1/(i + 1) + 1)
var pred = Fmt.f(8, t, 3)
System.print(" %(i+1) %(obs) %(pred)")
}
}
show.call(fib1000.call(), "First 1000 Fibonacci numbers:")</syntaxhighlight>
{{out}}
<pre>
First 1000 Fibonacci numbers:
Digit Observed Predicted
1 0.301 0.301
2 0.177 0.176
3 0.125 0.125
4 0.096 0.097
5 0.080 0.079
6 0.067 0.067
7 0.056 0.058
8 0.053 0.051
9 0.045 0.046
</pre>
=={{header|XPL0}}==
The program is run like this: benford < fibdig.txt
Luckly, there's no need to show how the first digits of the Fibonacci sequence were obtained.
<syntaxhighlight lang "XPL0">int D, N, Counts(10);
real A, E;
[for D:= 0 to 9 do Counts(D):= 0;
for N:= 1 to 1000 do
[D:= ChIn(1) & $0F; \ASCII to binary
Counts(D):= Counts(D)+1;
];
Text(0, "Digit Actual Expected Difference^m^j");
for D:= 1 to 9 do
[IntOut(0, D); ChOut(0, ^ );
A:= float(Counts(D))/1000.;
RlOut(0, A);
E:= Log(1. + 1./float(D));
RlOut(0, E);
RlOut(0, E-A);
CrLf(0);
];
]</syntaxhighlight>
{{out}}
<pre>
Digit Actual Expected Difference
1 0.30100 0.30103 0.00003
2 0.17700 0.17609 -0.00091
3 0.12500 0.12494 -0.00006
4 0.09600 0.09691 0.00091
5 0.08000 0.07918 -0.00082
6 0.06700 0.06695 -0.00005
7 0.05600 0.05799 0.00199
8 0.05300 0.05115 -0.00185
9 0.04500 0.04576 0.00076
</pre>
=={{header|Yabasic}}==
Using function from https://www.rosettacode.org/wiki/Fibonacci_sequence#Yabasic
<syntaxhighlight lang="vb">n = 1000
dim actual(n)
for nr = 1 to n
num$ = str$(fibonacci(nr))
j = val(left$(num$,1))
actual(j) = actual(j) + 1
next
print "First 1000 Fibonacci numbers"
print "Digit ", "Actual freq ", "Expected freq"
for i = 1 to 9
freq = frequency(i)*100
print i using("###"), " ", (actual(i)/10) using("##.###"), " ", freq using("##.###")
next
end
sub frequency(n)
return (log10(n+1) - log10(n))
end sub
sub log10(n)
return log(n) / log(10)
end sub
sub fibonacci(n)
local n1, n2, k, sum
n1 = 0
n2 = 1
for k = 1 to abs(n)
sum = n1 + n2
n1 = n2
n2 = sum
next k
if n < 0 then
return n1 * ((-1) ^ ((-n) + 1))
else
return n1
end if
end sub</syntaxhighlight>
{{out}}
<pre>First 1000 Fibonacci numbers
Digit Actual freq Expected freq
1 30.100 30.103
2 17.700 17.609
3 12.500 12.494
4 9.600 9.691
5 8.000 7.918
6 6.700 6.695
7 5.600 5.799
8 5.300 5.115
9 4.500 4.576</pre>
=={{header|zkl}}==
{{trans|Go}}
<
(0).pump(1000,List,fcn(ab){ab.append(ab.sum(0.0)).pop(0)}.fp(L(1,1))),
"First 1000 Fibonacci numbers");
Line 2,604 ⟶ 4,871:
(1.0+1.0/i).log10()))
}
}</
{{trans|CoffeeScript}}
<
fcn fibgen(a,b) { return(a,self.fcn.fp(b,a+b)) } //-->L(fib,fcn)
Line 2,626 ⟶ 4,893:
println("Leading digital distribution of the first 1,000 Fibonacci numbers");
println("Digit\tActual\tExpected");
foreach i in ([1..9]){ println("%d\t%.3f\t%.3f".fmt(i,actual[i], expected[i-1])); }</
{{out}}
<pre>
Line 2,641 ⟶ 4,908:
9 0.045 0.046
</pre>
=={{header|ZX Spectrum Basic}}==
{{trans|Liberty BASIC}}
<
20 DIM b(9)
30 LET n=100
Line 2,667 ⟶ 4,933:
1060 NEXT j
1070 RETURN
</syntaxhighlight>
The results obtained are adjusted fairly well, except for the number 8. This occurs with Sinclair BASIC, Sam BASIC and SpecBAS fits.
| |||