Monte Carlo simulation Lotto Germany

Question

Currently, I am working on programming a Monte Carlo simulation to investigate the likelihood of the rarest number (13) in German lottery draws being drawn only 361 times or less in 3570 draws.

import random as rnd
import numpy as np

#Variablen setzen 
anzahl_sim=10000
lottozahlen_alle = list(range(1,50))
sim_wie_realität = 0
liste_niedrgiste_häufigkeit=[]

#Schleife für Simulation
for i in range(anzahl_sim):
    häufigkeit_zahlen_aus_lottoziehung = {zahl: 0 for zahl in lottozahlen_alle}
    for j in range(3570):           #sci.py.stats.hypergeom für matrixbasierte Simulatione
        x = rnd.sample(lottozahlen_alle, 6)     #rnd.sample funktioniert wie "Urnenziehung"
        for zahl in x:
                häufigkeit_zahlen_aus_lottoziehung[zahl] += 1
    summe_häufigkeit_bei_3570_ziehungen = sum(häufigkeit_zahlen_aus_lottoziehung.values())
    niedrigste_zahl = min(häufigkeit_zahlen_aus_lottoziehung, key=häufigkeit_zahlen_aus_lottoziehung.get)
    niedrigste_häufigkeit = häufigkeit_zahlen_aus_lottoziehung[niedrigste_zahl]
    if niedrigste_häufigkeit <= 361:      #Prüfen, ob Simulation von 3570 Lottoziehungen Ergebnis wie Realität aufweist
        sim_wie_realität += 1
    liste_niedrgiste_häufigkeit.append(niedrigste_häufigkeit)   #Sammlung von seltensten Kugeln bei 3570 Ziehungen

array_niedrgiste_häufigkeit=np.array(liste_niedrgiste_häufigkeit)

p_sim_wie_realität = sim_wie_realität/anzahl_sim
avg_experiment = array_niedrgiste_häufigkeit.sum(axis=0)/anzahl_sim
std_experiment = array_niedrgiste_häufigkeit.std(axis=0)
var_experiment = array_niedrgiste_häufigkeit.var(axis=0)
entfernung_in_std = abs((361-avg_experiment)/std_experiment)

print("Mittelwert von Anzahl niedrigste Ziehung: ", avg_experiment)
print("Standardabweichung: ", std_experiment)
print("Varianz: ", var_experiment)
print("Entfernung in Standardabweichung realtität von Mittelwert Simulation: ", entfernung_in_std)
print("P(eine Zahl ist seltener als 361 mal oder 361 mal gezogen worden) für ", anzahl_sim, "Simmulationen: ", p_sim_wie_realität)

I have already thought a lot about the problem and believe that the following Monte Carlo simulation accurately examines this event's probability. However, I find it difficult to precisely calculate the theoretical probability beyond the Monte Carlo simulation. Perhaps you have a formula for me and know how I can optimize my code for the Monte Carlo simulation, especially since the random.sample function consumes a lot of computing time.