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('{:1d} {:6.4f} {:6.4f}'.format(r, p[r], pbin[r]))
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.