Benford's law: Difference between revisions

=={{header|11l}}==

{{trans|D}}

V a = 1.0

V b = 1.0

 

L(p) benford(get_fibs())

{{out}}

<pre>

 

=={{header|8th}}==

: n:log10e ` 1 10 ln / ` ;

 

benford-test

bye

 

{{out}}

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

Ada.Text_IO.New_Line;

end loop;

 

{{out}}

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.

 

=={{header|Aime}}==

sum(text a, text b)

{

 

0;

{{out}}

<pre> expected found

{{works with|ALGOL 68G|Any - tested with release 2.8.3.win32}}

Uses Algol 68G's LONG LONG INT which has programmer specifiable precision.

# set the number of digits for LONG LONG INT values #

PR precision 256 PR

# compare to the probabilities expected by Benford's law #

compare to benford( digit probability( fn ) )

{{out}}

<pre>

=={{header|AutoHotkey}}==

{{works with|AutoHotkey_L}}(AutoHotkey1.1+)

fib := NStepSequence(1, 1, 2, 1000)

Out := "Digit`tExpected`tObserved`tDeviation`n"

}

return, (Carry ? Carry : "") . Result

'''Output:'''

<pre>Digit Expected Observed Deviation

 

=={{header|AWK}}==

# syntax: GAWK -f BENFORDS_LAW.AWK

BEGIN {

function abs(x) { if (x >= 0) { return x } else { return -x } }

function log10(x) { return log(x)/log(10) }

{{out}}

<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)

GET "libhdr"

 

RESULTIS 0

}

{{Out}}

<pre>

 

=={{header|C}}==

#include <stdlib.h>

#include <math.h>

 

return EXIT_SUCCESS;

 

{{Out}}

 

=={{header|C++}}==

//if used prepackaged you can compile writing "g++ -std=c++11 -lcln yourprogram.cpp -o yourprogram"

#include <cln/integer.h>

return 0 ;

}

{{out}}

<pre> found expected

 

=={{header|Clojure}}==

(:gen-class))

 

(show-benford-stats y))

 

{{Output}}

<pre>

 

=={{header|CoffeeScript}}==

a = 1; b = 0

return () ->

console.log "Digit\tActual\tExpected"

for i in [1..9]

{{out}}

<pre>Leading digital distribution of the first 1,000 Fibonacci numbers

 

=={{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

(map 'list #'list '(1 2 3 4 5 6 7 8 9)

actual-distribution

 

<pre>; *fib1000* is a list containing the first 1000 numbers in the Fibonnaci sequence

=={{header|Crystal}}==

{{trans|Ruby}}

 

EXPECTED = (1..9).map{ |d| Math.log10(1 + 1.0 / d) }

# just to show that not all kind-of-random sets behave like that

 

{{out}}

=={{header|D}}==

{{trans|Scala}}

 

double[2][9] benford(R)(R seq) if (isForwardRange!R && !isInfinite!R) {

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

===Alternative Version===

The output is the same.

std.algorithm, std.array;

 

f * 100.0 / N, expected[i] * 100,

abs((f / double(N)) - expected[i]) * 100);

=={{header|Delphi}}==

See [https://rosettacode.org/wiki/Benford%27s_law#Pascal Pascal].

=={{header|Elixir}}==

def distribution(n), do: :math.log10( 1 + (1 / n) )

end

 

 

{{out}}

 

=={{header|Erlang}}==

-module( benfords_law ).

-export( [actual_distribution/1, distribution/1, task/0] ).

[Key | _] = erlang:integer_to_list( N ),

dict:update_counter( Key - 48, 1, Dict ).

{{out}}

<pre>

For Fibonacci code, see https://rosettacode.org/wiki/Fibonacci_sequence#F.89

 

 

let fibonacci = Seq.unfold (fun (x, y) -> Some(x, (y, x + y))) (0I,1I)

printfn "\nBenford's law for 1 through 9:"

benfordLawFigures |> List.iter (fun f -> printf $"{f:N5} ")

 

{{out}}

 

=={{header|Factor}}==

kernel math math.functions math.statistics math.text.utils

sequences ;

.header leading histogram [ .digit-report ] assoc-each ;

{{out}}

<pre>

 

=={{header|Forth}}==

: f2drop fdrop fdrop ;

 

set-precision ;

 

{{Out}}

<pre>Gforth 0.7.2, Copyright (C) 1995-2008 Free Software Foundation, Inc.

=={{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

 

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

 

 

! compute most significant digits, Fibonacci like.

implicit none

write(6,*) (benfordsLaw(i),i=1,9),'THE LAW'

write(6,*) (count(i)/1000.0 ,i=1,9),'LEADING FIBONACCI DIGIT'

 

=={{header|FreeBASIC}}==

{{libheader|GMP}}

' compile with: fbc -s console

 

Print : Print "hit any key to end program"

Sleep

{{out}}

<pre>First 1000 Fibonacci numbers

 

=={{header|Go}}==

 

import (

math.Log10(1+1/float64(i+1)))

}

{{Out}}

<pre>

Uses [[Fibonacci_sequence#Analytic_8|Fibonacci sequence analytic formula]]

{{trans|Java}}

def population = (0..<size).collect { generator(it) }

def firstDigits = [0]*10

}

firstDigits

 

'''Test:'''

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))

 

'''Output:'''

 

=={{header|Haskell}}==

import Data.Char (digitToInt)

 

| d <- [1 .. 9] ]

 

{{out}}

<pre>[(1,0.301,0.301029995663981),

The following solution works in both languages.

 

 

procedure main()

if n%2 = 1 then return [c+d, d]

else return [d, c]

 

Sample run:

=={{header|J}}==

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

assert '0.9999' -: 6j4 ": (normalize TALLY_BY_KEY) r benford >: i.9

 

 

=={{header|Java}}==

import java.util.Locale;

 

}

}

The output is:

<pre>1 0.301000 0.301030

 

=={{header|JavaScript}}==

.reduce(

(fib, _, i) => i < 2 ? (

]);

 

{{output}}

<pre>0: (3) [1, 0.301, 0.3010299956639812]

This implementation shows the observed and expected number of occurrences together with the χ² statistic.

 

# Numerical accuracy is insufficient beyond about 1450.

def fibonacci(n):

;

 

{{out}}

<pre>First 100 fibonacci numbers:

=={{header|Julia}}==

 

 

ben(l) = [count(x->x==i, map(n->string(n)[1],l)) for i='1':'9']./length(l)

 

{{Out}}

<pre>julia> benford([fib(big(n)) for n = 1:1000])

 

=={{header|Kotlin}}==

 

interface NumberGenerator {

}

 

 

=={{header|Liberty BASIC}}==

Using function from

http://rosettacode.org/wiki/Fibonacci_sequence#Liberty_BASIC

dim bin(9)

 

fiboI = a

end function

 

{{out}}

 

=={{header|Lua}}==

expected = {}

for i = 1, 9 do

for i = 1, 9 do

print(i, actual[i] / n, expected[i])

{{out}}

<pre>digit actual expected

 

=={{header|Mathematica}} / {{header|Wolfram Language}}==

Table[{d, N@Count[fibdata, d]/Length@fibdata, Log10[1. + 1/d]}, {d, 1,

{{out}}

<pre>1 0.301 0.30103

 

=={{header|NetRexx}}==

options replace format comments java crossref symbols nobinary

 

brenfordDeveation(fibList)

return

{{out}}

<pre>

 

=={{header|Nim}}==

import strformat

 

let distrib = actualDistrib(fibSeq)

echo "Fibonacci numbers first digit distribution:\n"

 

{{out}}

=={{header|Oberon-2}}==

{{Works with|oo2c version 2}}

MODULE BenfordLaw;

IMPORT

END

END BenfordLaw.

{{out}}

<pre>

https://rosettacode.org/wiki/Fibonacci_sequence#Arbitrary_Precision<br>

Note the remark about the compilation of the program there.

open Num

 

List.iter (Printf.printf "%f ") (List.map benfords_law xvalues) ;

Printf.printf "\n" ;;

{{out}}

<pre>

 

=={{header|PARI/GP}}==

my(t=vector(9,n,sum(i=1,#v,v[i]==n)));

print("Digit\tActual\tExpected");

lucas(n)=fibonacci(n-1)+fibonacci(n+1);

dist(fibonacci)

{{out}}

<pre>Digit Actual Expected

 

=={{header|Pascal}}==

{$IFDEF FPC}{$MODE Delphi}{$ELSE}{$APPTYPE CONSOLE}{$ENDIF}

uses

writeln(i:5,dgtCnt[i]:7,expectedCnt[i]:10,reldiff:10:5,' %');

end;

 

<pre>Digit Count Expected rel Diff

 

=={{header|Perl}}==

use strict ;

use warnings ;

$result = sprintf ( "%.2f" , 100 * $expected[ $i - 1 ] ) ;

printf "%15s %%\n" , $result ;

{{Out}}

<pre>

=={{header|Phix}}==

{{trans|Go}}

<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: #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>

{{out}}

<pre>

 

=={{header|Picat}}==

 

N = 1000,

foreach(E in List)

Map.put(E, cond(Map.has_key(E),Map.get(E)+1,1))

 

{{out}}

=={{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.

(scl 4)

(load "@lib/misc.l")

 

(bye)

{{Out}}

<pre>

 

=={{header|PL/I}}==

(fofl, size, subrg):

Benford: procedure options(main); /* 20 October 2013 */

end tally;

end Benford;

Results:

<pre>

 

=={{header|PL/pgSQL}}==

WITH recursive

constant(val) AS

(select cast(corr(probability_theoretical,probability_real) as numeric(5,4)) correlation

from benford) c

 

=={{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.

$url = "https://oeis.org/A000045/b000045.txt"

$file = "$env:TEMP\FibonacciNumbers.txt"

 

Remove-Item -Path $file -Force -ErrorAction SilentlyContinue

{{Out}}

<pre>

{{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?

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

findall(B, (between(1,9,N), benford(N,B)), Benford),

findall(C, firstchar(C), Fc), freq(Fc, Freq),

{{out}}

<pre>?- go.

 

=={{header|PureBasic}}==

NewMap d1.i()

Dim fi.s(#MAX_N)

 

PrintN(~"\nPress Enter...")

{{out}}

<pre>Dig. Cnt. Exp. Dif.

=={{header|Python}}==

Works with Python 3.X & 2.7

from itertools import islice, count

from collections import Counter

 

# just to show that not all kind-of-random sets behave like that

{{out}}

<pre>fibbed Benfords deviation

 

=={{header|R}}==

pbenford <- function(d){

return(log10(1+(1/d)))

 

print(data)

{{out}}

digit obs.frequency exp.frequency dev_percentage

 

=={{header|Racket}}==

 

(define (log10 n) (/ (log n) (log 10)))

;; 7: 5.8% 5.8%

;; 8: 5.1% 5.1%

 

=={{header|Raku}}==

(formerly Perl 6)

{{Works with|rakudo|2016-10-24}}

 

sub show(%distribution) {

 

multi MAIN($file) { show benford $file.IO.lines }

 

'''Output:''' First 1000 Fibonaccis

For the extra credit stuff, it was chosen to generate Fibonacci and factorials rather than find a

web─page with them listed, &nbsp; as each list is very easy to generate.

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. */

do f=1 for 9; say pad center(f,7) pad center(format(!.f/N,,length(N-2)),w1) #.f

end /*k*/

{{out|output|text=&nbsp; when using the default (1000 numbers) &nbsp; for the input:}}

<pre>

 

=={{header|Ring}}==

# Project : Benford's law

 

if y = 1 return 1 ok

if y > 1 return fibonacci(y-1) + fibonacci(y-2) ok

Output:

<pre>

=={{header|Ruby}}==

{{trans|Python}}

 

def fib(n)

 

# just to show that not all kind-of-random sets behave like that

 

{{out}}

 

=={{header|Run BASIC}}==

N = 1000

for i = 0 to N - 1

next i

end function

<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>

 

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.

 

extern crate num_traits;

extern crate num;

 

}

{{out}}

<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 }

 

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>

=={{header|Scheme}}==

{{works with|Chez Scheme}}

 

(define benford-probability

(display-table "Rnd/1T/1M" (lambda () (1+ (random 1000000000000))) 1000000)

(let ((craters (list->vector (list-read-file "moon_craters.lst"))))

{{out}}

The first one thousand Fibonnaci numbers.

 

=={{header|Sidef}}==

var fibonacci = 1000.of {|i| fib(i).digit(0) }

 

"%.2f".sprintf(100 * expected[i - 1]),

)

 

{{out}}

The query is the same for any number sequence you care to put in the <tt>benford</tt> table.

 

create table benford (num integer);

 

 

-- Tidy up

 

{{out}}

=={{header|Stata}}==

 

set obs 1000

scalar phi=(1+sqrt(5))/2

8 | 53 51.15252245 |

9 | 45 45.75749056 |

 

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]:

 

chisq=sum((f-p):^2:/p)

chisq

chi2tail(8,chisq)

.9999942179

 

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:

 

set obs 500

scalar s=2+sqrt(2)

chi2tail(8,chisq)

.0387287805

 

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}}==

 

 

/* Reads from a file and returns the content as a String */

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))

{{out}}

<pre>$ ./Benford

 

=={{header|Tcl}}==

# Count the leading digits (RE matches first digit in each number,

# even if negative)

[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]}]]

return $result

}

{{out}}

<pre>

 

=={{header|VBA (Visual Basic for Application)}}==

Sub BenfordLaw()

 

 

}

 

=={{header|Visual FoxPro}}==

#DEFINE CTAB CHR(9)

#DEFINE COMMA ","

RETURN LOG10(1 + 1/d)

ENDFUNC

{{out}}

<pre>

{{trans|Go}}

{{libheader|Wren-fmt}}

 

fn fib1000() []f64 {

println(" ${i+1} ${f64(n)/f64(c.len):9.3f} ${math.log10(1+1/f64(i+1)):8.3f}")

}

 

{{out}}

{{trans|Go}}

{{libheader|Wren-fmt}}

 

var fib1000 = Fn.new {

}

 

 

{{out}}

=={{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");

(1.0+1.0/i).log10()))

}

{{trans|CoffeeScript}}

 

fcn fibgen(a,b) { return(a,self.fcn.fp(b,a+b)) } //-->L(fib,fcn)

println("Leading digital distribution of the first 1,000 Fibonacci numbers");

println("Digit\tActual\tExpected");

{{out}}

<pre>

=={{header|ZX Spectrum Basic}}==

{{trans|Liberty BASIC}}

20 DIM b(9)

30 LET n=100

1060 NEXT j

1070 RETURN

The results obtained are adjusted fairly well, except for the number 8. This occurs with Sinclair BASIC, Sam BASIC and SpecBAS fits.