Příklad využití Monte Carlo simulace pro testování statistické významnosti korelace

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stat

rng=np.random.default_rng()
print(rng)

Generator(PCG64)

# načtení dat
data=np.loadtxt('data.txt')
x=data[:,0] #vyska
y=data[:,1] #vaha

# graf hodnot
fig,ax=plt.subplots()
ax.scatter(x,y)
ax.set_xlabel('x',fontsize=12)
ax.set_ylabel('y',fontsize=12)
plt.show()

#Pearsonův korelační koeficient
Pearson,P=stat.pearsonr(x,y)
print('korelace = ',Pearson)

korelace =  0.6513609129813575

#shuffling
N=10000 #počet simulací
i_counter=0
rhos=np.empty(N)
for i in range(N):
    xs=rng.permutation(x)
    rhos[i]=np.corrcoef(xs,y)[0,1]
    if np.abs(rhos[i])>=np.abs(Pearson):
        i_counter+=1

print('P value (calculated) = ',P)
print('P value (from simulation) = ',i_counter/N)

P value (calculated) =  0.00012968327111122398
P value (from simulation) =  0.0004

# histoghram korelací
fig,ax=plt.subplots()
ax.hist(rhos,bins=10)
ax.plot([Pearson,Pearson],[0,2500],c='red')
ax.set_xlabel('korelace',fontsize=12)
ax.set_ylabel('počet případů',fontsize=12)
plt.show()