Using a for loop to count the digits of pi

If the digits of π were random, then we would expect each of the integers to occur with approximately equal frequency in the decimal representation of π. In this problem we are going to see if the digits of π do indeed appear to be random.

The first line of your script stores the first 100 digits of π in the vector pi_digits.

1. Create a vector f50 such that f50(i) is equal to the frequency with which the integer i-1 appears among the first 50 digits of π.
2. Replicate it for f100.

I have the following code:

pi_digits=[3 1 4 1 5 9 2 6 5 3 5 8 9 7 9 3 2 3 8 4 6 2 6 4 3 3 8 3 2 7 9 5 0 2 8 8 4 1 9 7 1 6 9 3 9 9 3 7 5 1 0 5 8 2 0 9 7 4 9 4 4 5 9 2 3 0 7 8 1 6 4 0 6 2 8 6 2 0 8 9 9 8 6 2 8 0 3 4 8 2 5 3 4 2 1 1 7 0 6 7];
f50=zeros(1,10);
for i=0:9
    f50(i+1)=sum(pi_digits(1:50)==i)
end
f100=zeros(1,10);
for i=0:9
    f100(i+1)=sum(pi_digits(1:100)==i)
end

the output returns the correct frequency vector, yet I get an error that it is the incorrect value.

[Walter Roberson]

I recommend that you semi-colon to the end of the sum() lines, as you do not need the intermediate values displayed.

[Martin]

Do you have any suggestions as to why I receive an incorrect value when the output is correctly counting the frequency of integers 0-9?

[檮杌]

Which part do you think is incorrect? You can use histcounts to count frequencies of vector entities.

pi_digits=[3 1 4 1 5 9 2 6 5 3 5 8 9 7 9 3 2 3 8 4 6 2 6 4 3 3 8 3 2 7 9 5 0 2 8 8 4 1 9 7 1 6 9 3 9 9 3 7 5 1 0 5 8 2 0 9 7 4 9 4 4 5 9 2 3 0 7 8 1 6 4 0 6 2 8 6 2 0 8 9 9 8 6 2 8 0 3 4 8 2 5 3 4 2 1 1 7 0 6 7];
f50=zeros(1,10);
for i=0:9
    f50(i+1)=sum(pi_digits(1:50)==i);
end
f100=zeros(1,10);

for i=0:9
    f100(i+1)=sum(pi_digits(1:100)==i);
end

histcount50 = histcounts(pi_digits(1:50))

histcount50 = 1×10

1 5 5 9 4 5 4 4 5 8

histcount100 = histcounts(pi_digits)

histcount100 = 1×10

8 8 12 12 10 8 9 8 12 13

isequal(f50, histcount50)

ans = logical
   1

isequal(f100, histcount100)

ans = logical
   1

[Sam Chak]

@Martin (OP): "In this problem we are going to see if the digits of π do indeed appear to be random."

Sounds interesting. I believe this is still an open problem in Mathematics and Algorithmics. I heard many of my math friends say its decimal expansion goes on forever without repeating, but no one has actually tested it 'forever'.

Do you think the first 100 digits of π are sufficient to demonstrate the randomness?

[Martin]

@檮杌 I believe my code is correct but it is incorrect on matlab grader.

[Martin]

@Sam Chak While I do not believe the first 100 digitts are sufficient to desmonstrate the randomness, I do want to get the marks! Yet while the frequency vectors are correct as shown by the histcount, it still says it is incorrect.

[Sam Chak]

No worries @Martin, I'm just mathematically curious! It will take some time to troubleshoot the issue. Could that be a technical problem on MATLAB Grader itself?

[William Rose]

@Martin,

[edit: I attached the text file of a million decimal digits (plus the initial "3.") from here, in case you are having trouble downloading it.]

if you're interested, you can download pi to a million decimal digits here. The site already reports the histogram of digits, so you can check your algorithm. Have fun!

[Walter Roberson]

@William Rose

I was not able to follow the ftp link, at least not with Firefox. I will try with a command line tool.

[Martin]

@Sam Chak The question itself is definitely interesting, but yeah that's what I'm thinking as well... to office hours I go...

[Walter Roberson]

Nope, cannot resolve the site name wuarchive.wustl.edu

[Sam Chak]

I think wuarchive-ftpd was defunct.

https://en.wikipedia.org/wiki/WU-FTPD

Can find some archives here:

http://ftp.ist.utl.pt/pub/ftp/servers/wuarchive-ftpd/wu-ftpd-attic/wustl.edu/

[Sam Chak]

Can also find something interesing here:

https://pi2e.ch/blog/2017/03/10/pi-digits-download/

Up to 22.4 Trillion Digits of π!

[Dyuman Joshi]

@Sam, but the discussion here is for pi, not the factorial of pi :P

[Sam Chak]

@Dyuman Joshi, I am "irrational", pun intended.

Answer 1

[Daniel]

MA207?

need to divide f50 by 50 and f100 by 100.

apparently it wants the answer to be in decimal.

[Walter Roberson]

Ah, the original Question does ask about "frequency" rather than about counts, so I can see why they might normalize by the number of entries.

[Sam Chak]

Thank you, @Daniel. I interpreted "frequency" as the number of occurrences within a given 100-digit π number. I totally overlooked that the "frequency" in question is somehow mathematically defined in statistics.

https://en.wikipedia.org/wiki/Frequency_(statistics)

Answer 2

[Sam Chak]

Hi @Martin

Please try comparing the histogram with yours.

num2str(pi, 1000);
digits(100);    % show 100 digits of π
numPi = vpa(pi)

numPi = 

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

c     = char(numPi);

% Find the decimal point:
pos   = strfind(c, '.')

pos = 2

% Begin to count after the decimal point
d = arrayfun(@str2num, c(pos+1:end));

% Plot the histogram for comparing with your Answer in MATLAB Grader
histogram(d, 10);
xt = linspace(0.5, 8.6, 10);
xticks(xt);
xticklabels({'0', '1', '2', '3', '4', '5', '6', '7', '8', '9'})
xlabel('Decimal Digit')
ylabel('Frequency')
title('Distribution of the digits of \pi');

[Sam Chak]

Hi @Martin

This code converts the number π into an array of integers and counts the frequency of the digits. The code also verifies the results obtained by you. Please review the question to clarify whether you want to count the frequency of the digits in the array of integers or only the decimal digits. Additionally, if I use digits(100), the last digit will be rounded up.

num2str(pi, 1000);
digits(102);    % show 102 digits of π
numPi = vpa(pi)

numPi = 

3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798

c     = char(numPi);
c(2)  = []      % remove decimal point

c = '314159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798'

% Count the digit from array of integers: 3, 1, 4, ...
d = arrayfun(@str2num, c(1:end-2))  % exclude the 101st and 102nd position in the sequence

d = 1×100

3 1 4 1 5 9 2 6 5 3 5 8 9 7 9 3 2 3 8 4 6 2 6 4 3 3 8 3 2 7

histcount100 = histcounts(d)

histcount100 = 1×10

8 8 12 12 10 8 9 8 12 13

% Plot the histogram
histogram(d, 10);
xt = linspace(0.5, 8.5, 10);
xticks(xt);
xticklabels({'0', '1', '2', '3', '4', '5', '6', '7', '8', '9'})
xlabel('Digit')
ylabel('Frequency')
title('Distribution of the digits of \pi');

[Martin]

Thanks I appreciate it! After further review I am convinced there is an error in the solution on Matlab Grader.

[Stephen23]

Another approach using 1000 digits from https://pi2e.ch/blog/2017/03/10/pi-digits-download/

T = '31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989';
histogram(T(1:100)-'0',0:10)

Confirming a few random digits:

nnz(T(1:100)=='0')

ans = 8

nnz(T(1:100)=='6')

ans = 9

nnz(T(1:100)=='9')

ans = 13