Simulace odhadu pravděpodobnosti, že padne orel při házení mincí

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import comb  #kombinační číslo

N=10 #počet hodů

posteriorní hustota pravděpodobnosti $$f(p|k)= (N+1)\frac{N!}{(N-k)! k!}p^k (1-p)^{N-k}$$

def f(k,N,p):
    return (N+1)*comb(N,k)*p**k*(1-p)**(N-k)

simulace

p_plot=np.linspace(0,1,100) #x-ová souřadnice pro vykreslení grafu 
k=0 #počet orlů
for i in range(N):
    x=np.random.choice(2,1) #simulace jednoho hodu: 0...panna, 1..orel
    if x>0:  #padl orel
        k+=1
    fig,ax=plt.subplots()
    ax.plot(p_plot,f(k,i+1,p_plot),label='N = '+str(i+1)+', k = '+str(k))
    ax.legend()
    ax.set_xlabel('p')
    ax.set_ylabel('f(p|D)')
    ax.set_xlim(0,1)