https://math.stackexchange.com/q/3947387/650339
\(\pi\) is waarschijnlijk een normaal getal, betekenend dat de decimalen progressie is random is. De decimalen [0,1,2,3,4,5,6,7,8,9] komen allen even waarschijnlijk voor. Dit is getest voor vele miljoenen decimalen en is waarschijnlijk waar.
Onderstaand een methode om \(\pi\) te bepalen uit random getallen. Dus in principe hoe \(\pi\) te berekenen uit de decimalen van \(\pi\)!
We gaan uit van de PDF (probability density functie) van een normaal verdeling:
$$f(x)={\frac {1}{\sigma {\sqrt {2\pi }}}}e^{-{\frac {1}{2}}\left({\frac {x-\bar{x}}{\sigma }}\right)^{2}}$$
De dichtheid bij de gemiddelde waarde: \(\bar{x}\) is gelijk aan de voorkant van de PDF. Deze gemiddelde dichtheid kan worden benaderd door een discrete verdeling. Waar \(\bar{n}\) het aantal getelde / waargenomen waarden is bij de gemiddelde \(\bar{x}\), \(N\) is de steekproefomvang en \(\Delta x\) de afstand tussen twee discrete intervallen.
$$f(\bar{x})=\frac {1}{\sigma {\sqrt {2\pi }}} \approx \frac{\bar{n}}{\Delta x N}$$
Dus als we normaal verdeeld zijn, kunnen we \(\pi\) berekenen met:
$$\boxed{\pi \approx \frac{1}{2} \cdot \left( \frac{\Delta x N}{ \bar{n} \sigma} \right)^{2}}$$
Laten we een simulatie maken, we kiezen een cijfer uit [0,1,2,3,4,5,6,7,8,9] we nemen monsters "met vervanging". De standaarddeviatie van het gemiddelde kan worden berekend als: a = 0, b = 9 en s steekproefomvang. Een discrete uniforme verdeling:
$$\sigma=\sqrt{\frac{(b-a+1)^{2}-1}{12s}}$$
De afstand tussen discrete intervallen kan worden berekend met (empirisch bevonden):
$$\Delta x=\frac{1}{s}$$
Nu zijn we klaar om bijvoorbeeld een simulatie te doen: neem s = 800 willekeurige steekproeven van een van: [0,1,2,3,4,5,6,7,8,9] ("met vervanging") en bepaal de gemiddelde waarde. We herhalen elk monster met N = 150000 proeven, dit is het aantal steekproefgemiddelden. Van een enkele proef kunnen we tellen hoeveel monsters de verwachte gemiddelde waarde hebben: \(\bar{x}=4.5\) (alleen als de steekproefomvang even is). Hieruit kunnen we de geschatte \(\pi\) in één keer berekenen. We herhalen elke proef 10.000 keer, zodat we het gemiddelde van \(\pi\) kunnen bepalen.
Linker grafiek: Alle herhaalde pogingen, de grijze stippen zijn de willekeurig gegenereerde waarden. De blauwe lijn is de theoretische verdeling.
Middelste grafiek: de gemiddelde getelde waarnemingen n bij de verwachte waarde: \(\bar{x}=4.5\). Dit geeft een normale verdeling en de standaarddeviatie van de gemiddelde telling kan worden bepaald.
Rechter grafiek: de waarde van \(\pi\) kan bepaald door de omkaderde formule, merk op dat de resulterende verdeling niet normaal is. De boven- en ondergrens 95% wordt bepaald uit de middelste grafiek.
Code: Selecteer alles
import numpy as np
import pandas as pd
#Samples, trials (pick trials >1000 to avoid zero counts) and repeats
samples=100
trials=2000
repeats=2000
#Set arrays to zero
piarray=np.zeros(repeats)
narray=np.zeros(repeats)
#Detemermine: Stdev, And dx
var=((10-1+1)**2-1)/12
stdev=np.sqrt(var)
stdevt=stdev/np.sqrt(samples)
dx=1/samples
for p in range(repeats):
#Create Array random numbers
random=np.random.choice([0,1,2,3,4,5,6,7,8,9],[trials,samples])
#Determine mean from random no of numbers (of samples)
m=np.mean(random,axis=1)
mn=np.mean(m)
#Create dataframe and count number of observation per interval
df=pd.DataFrame({'m' : m})
dfg=df.groupby(['m'])['m'].agg(['count']).reset_index()
#Histogram x=bins and y is density
x=dfg['m'].to_numpy()
y=dfg['count'].to_numpy()
y=y/(dx*trials)
#Count number of observation on mean: 4.5
out=dfg[(dfg['m'] < 4.5 +dx/2) & (dfg['m'] > 4.5-dx/2)]
nc=out['count'].to_numpy()
n=np.sum(nc)
narray[p]=n
#Calculate pi and fill array
pi=0.5*(trials*dx/(stdevt*n))**2
piarray[p]=pi
nmean=np.mean(narray)
nstdev=np.std(narray)
nstdevm=nstdev/np.sqrt(repeats)
pim=0.5*(dx*trials/(stdevt*nmean))**2
piminc=0.5*(dx*trials/(stdevt*(nmean+2*nstdevm)))**2
pimaxc=0.5*(dx*trials/(stdevt*(nmean-2*nstdevm)))**2
print('pi: ' + str(np.round(pim,4)))
print(str(np.round(piminc,4)) + '<pi<' + str(np.round(pimaxc,4)) + ' (95%)')Ik merkte later dat \(\pi\) ook kan worden berekend met de dichtheid op de positie van \(\sigma\) (of een andere positie). Maar dan heb je het nummer \(e\) van Euler nodig. Dus in principe kan na het bepalen van \(\pi\) het getal van Euler bepaald worden.
$$\boxed{\pi \approx \frac{1}{2e} \cdot \left( \frac{\Delta x N}{ n_{\sigma} \sigma} \right)^{2}}$$
