Function Calculates the Probability of a Fibonacci Number

My task is to write a function that calculates the probability the the k-th digit of a Fibonacci Number in B is d. The elements of Prob are the probabilities of the k-th digit of a Fibonacci Number in B being d, for d = 0, 1, 2, ..., 9. n and k are positive integers and Prob is a 1x10 vector of probabilities. A is the set of the first n Fibonacci Numbers and B is the subset of the Fibonacci Numbers in A that have at least k digits. d is each of the digits 0-9. I have a small program written up, however I am stuck on where to go next. I know Prob has to be incoorporated somewhere, but I don't know where. Here is what I have so far:

function Prob = ex(n,k)
%
%
Prob = zeros(1,10);
A = zeros(1,n);
B = zeros(1,k);
position = 1;

while position <= n
    if n < 3
        Fn = 1;
    else
        Fn_2 = 1;
        Fn_1 = 1;
        for k = 3:n
            Fn = Fn_1 + Fn_2;
            Fn_2 = Fn_1;
            Fn_1 = Fn;
        end
    end
    position = position + 1;
end

% need to input d somewhere 

Answer 1

[Satyam]

Hi Sophie,

I understand that you're working on calculating the probability of the k-th digit of a Fibonacci number in a subset being a specific digit. Here is a step by step method to assist you in calculating the probability:

• Generate the first 'n' Fibonacci numbers accurately. Use a loop to iterate and calculate each Fibonacci number, storing them in an array 'A'.
• Create a subset 'B' of these numbers that have at least 'k' digits.
• For each number in your filtered subset, check the k-th digit and count its occurrences for each digit (0-9). Store these counts in the probability array 'Prob'.
• Convert these counts into probabilities by dividing by the total number of numbers in your subset 'B'.

Here is a code snippet modelling this approach:

function Prob = ex(n, k)
    % Initialize probability vector
    Prob = zeros(1, 10);
    
    % Generate first n Fibonacci numbers
    A = zeros(1, n);
    A(1) = 1;
    A(2) = 1;
    for i = 3:n
        A(i) = A(i-1) + A(i-2);
    end
    
    % Filter Fibonacci numbers with at least k digits
    B = [];
    for i = 1:n
        if length(num2str(A(i))) >= k
            B = [B, A(i)];
        end
    end
    
    % Calculate probabilities
    for i = 1:length(B)
        numStr = num2str(B(i));
        if length(numStr) >= k
            digit = str2double(numStr(k));
            Prob(digit + 1) = Prob(digit + 1) + 1;
        end
    end
    
    % Normalize to get probabilities
    if ~isempty(B)
        Prob = Prob / length(B);
    end
    disp(Prob)
end

ex(100,3);

    0.0899    0.1348    0.1011    0.1124    0.0674    0.1011    0.0899    0.1573    0.0562    0.0899

I hope this addresses your issue.

[John D'Errico]

@Satyam

Note that while this will work to an extent, it will fail, even for the example you show.

You have done all of these computations in double precision. But the 78th Fibonacci number is the largest one that is exactly representable as a double.

format long g
fibonacci([77 78 79])

ans = 1×3

1.0e+00 * 5.52793970088476e+15 8.94439432379146e+15 1.44723340246762e+16

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

flintmax('double')

ans = 
      9.00719925474099e+15

flintmax returns the largest number that is exactly representable as an integer for that precision. In the case of a double, it is 2^53-1. So while a double can encode numbers larger than that value, the low order digits of numbers beyond that point will in general not be correct, just random floating point trash.

And that means the lower order digits you will compute with this function are complete and utter garbage, beyond the 78'th Fibonacci number.

[Walter Roberson]

function Prob = ex(n, k)
    % Initialize probability vector
    Prob = zeros(1, 10);
    
    % Generate first n Fibonacci numbers
    A = zeros(1, n, 'sym');
    A(1) = 1;
    A(2) = 1;
    for i = 3:n
        A(i) = A(i-1) + A(i-2);
    end
    
    % Filter Fibonacci numbers with at least k digits
    B = [];
    for i = 1:n
        if length(char(A(i))) >= k
            B = [B, A(i)];
        end
    end
    
    % Calculate probabilities
    for i = 1:length(B)
        numStr = char(B(i));
        if length(numStr) >= k
            digit = str2double(numStr(k));
            Prob(digit + 1) = Prob(digit + 1) + 1;
        end
    end
    
    % Normalize to get probabilities
    if ~isempty(B)
        Prob = Prob / length(B);
    end
    disp(Prob)
end

ex(100,3);

    0.1236    0.0899    0.0899    0.0899    0.1011    0.0787    0.1011    0.1461    0.0787    0.1011