There are two stable isotopes of carbon, $\mathrm{^{12}C}$ and $\mathrm{^{13}C}$ (the radioactive $\mathrm{^{14}C}$ nucleus is present in nature in only trace amounts of the order of parts per trillion). Taking the abundance of $\mathrm{^{13}C}$ to be $x=0.0107$ (i.e. about 1%), we will calculate the relative amounts of buckminster fullerene, $\mathrm{C_{60}}$, with exactly zero, one, two, three and four $\mathrm{^{13}C}$ atoms. (This is important in nuclear magnetic resonance studies of fullerenes, for example, because only the $\mathrm{^{13}C}$ nucleus is magnetic and so detectable by NMR).

The number of $\mathrm{^{13}C}$ atoms in a population of carbon atoms sampled at random from a population with natural isotopic abundance follows a binomial distribution: the probability that, out of $n$ atoms, $m$ will be $\mathrm{^{13}C}$ (and therefore $n-m$ will be $\mathrm{^{12}C}$) is $$ p_m(n) = \binom{n}{m} x^m (1-x)^{n-m}. $$ We can, of course, calculate $p_m(60)$ exactly from this formula for $0 \le m \le 4$, but we can also simulate the sampling with the np.random.binomial method:

import numpy as np

n, x = 60, 0.0107
mmax = 4
m = np.arange(mmax + 1)

# Estimate the abundances by random sampling from the binomial distribution
ntrials = 10000
pbin = np.empty(mmax + 1)
for r in m:
    pbin[r] = np.sum(np.random.binomial(n, x, ntrials) == r) / ntrials

# Calculate and store the binomial coefficients nCm
nCm = np.empty(mmax + 1)
nCm[0] = 1
for r in m[1:]:
    nCm[r] = nCm[r - 1] * (n - r + 1) / r
# The "exact" answer from binomial distribution
p = nCm * x**m * (1 - x) ** (n - m)

print("Abundances of C60 as (13C)[m](12C)[60-m]")
print('m   "Exact"    Estimated')
print("-" * 24)
for r in m:
    print(f"{r:1d}   {p[r]:6.4f}      {pbin[r]:6.4f}")

For each value of r in the array m, we sample a large number of times (ntrial) from the binomial distribution described by $n=60$ and probability, $x=0.0107$. The comparison of these sample values with a given value of r yields a boolean array which can be summed (remembering that True evaluates to 1 and False evaluates to 0); division by ntrials then gives an estimate of the probability of exactly r atoms being of type $\mathrm{^{13}C}$ and the remainder of type $\mathrm{^{12}C}$.

The explicit loop over m could be removed by creating an array of shape (ntrials, mmax+1) containing all the samples, and summing over the first axis of this array in the comparison with the m array:

samples = np.random.binomial(n, x, (ntrials, mmax+1))
pbin = np.sum(samples == m, axis=0) / ntrials

The abundances of $\mathrm{^{13}C_m^{12}C_{60-m}}$ produced by our program are given as the following output.

Abundances of C60 as (13C)[m](12C)[60-m]
m   "Exact"    Estimated
------------------------
0   0.5244      0.5199
1   0.3403      0.3348
2   0.1086      0.1093
3   0.0227      0.0231
4   0.0035      0.0031

That is, almost 48% of $\mathrm{C_{60}}$ molecules contain at least one magnetic nucleus.