57 lines
1.4 KiB
Matlab
57 lines
1.4 KiB
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
|