In [63]:

    

%%latex
We define the binomial distribution $f$ as
\begin{align}
\Theta &= 0.011 \\
f(k;n,\Theta) &= {n \choose k}\Theta^k(1-\Theta)^{n-k}\\
\end{align}
Let $X$ be the measured isotope fractions, and $Y$ the real isotope fractions (i.e. if the natural abundance was zero).
\begin{align}
X_0 &= f(0;N,\Theta) \cdot Y_0 \\
X_1 &= f(1;N,\Theta) \cdot Y_0 + f(0;N-1,\Theta) \cdot Y_1 \\
X_2 &= f(2;N,\Theta) \cdot Y_0 + f(1;N-1,\Theta) \cdot Y_1 + f(0;N-2,\Theta) \cdot Y_2 \\
X_i &= \sum_{j=0}^{i} f(i-j;N-j,\Theta) \cdot Y_j
\end{align}

So we can define a matrix $F$, such that:
$$\forall i \ge j ~~~ F_{i,j} = f(i-j;N-j,\Theta) = {N-j \choose i-j} \Theta^{i-j} (1-\Theta)^{N-i}$$
and the rest are zeros. Then it is easy to see that $X = FY$ and since $F$ is non-singular we can solve $Y = F^{-1}X$.


    

    






We define the binomial distribution $f$ as
\begin{align}
\Theta &= 0.011 \\
f(k;n,\Theta) &= {n \choose k}\Theta^k(1-\Theta)^{n-k}\\
\end{align}
Let $X$ be the measured isotope fractions, and $Y$ the real isotope fractions (i.e. if the natural abundance was zero).
\begin{align}
X_0 &= f(0;N,\Theta) \cdot Y_0 \\
X_1 &= f(1;N,\Theta) \cdot Y_0 + f(0;N-1,\Theta) \cdot Y_1 \\
X_2 &= f(2;N,\Theta) \cdot Y_0 + f(1;N-1,\Theta) \cdot Y_1 + f(0;N-2,\Theta) \cdot Y_2 \\
X_i &= \sum_{j=0}^{i} f(i-j;N-j,\Theta) \cdot Y_j
\end{align}

So we can define a matrix $F$, such that:
$$\forall i \ge j ~~~ F_{i,j} = f(i-j;N-j,\Theta) = {N-j \choose i-j} \Theta^{i-j} (1-\Theta)^{N-i}$$
and the rest are zeros. Then it is easy to see that $X = FY$ and since $F$ is non-singular we can solve $Y = F^{-1}X$.

In [127]:

    

#######################################################
# Here you place the input data - counts for each mass 
# starting with M+0 (of course)
#
counts = [900, 100, 5, 900, 5000]
#
#######################################################

import numpy as np
from scipy.misc import comb
theta = 0.011

N = len(counts)
X = np.matrix(counts, dtype=float).T
X = X / X.sum()

F = np.matrix(np.zeros((N, N)))
for i in xrange(N):
    for j in xrange(i+1):
        F[i,j] = comb(N-j, i-j) * theta**(i-j) * (1-theta)**(N-j)
        
Y = F.I * X

print "The corrected isotope relative abundances are:"
print '-'*50
print '  |  '.join(map(lambda d: ' M + %d' % d, range(N)))
print '  |  '.join(map(lambda s: '%5.2f%%' % (s*100), Y.flat))
print '-'*50


    

    




The corrected isotope relative abundances are:
--------------------------------------------------
 M + 0  |   M + 1  |   M + 2  |   M + 3  |   M + 4
13.78%  |   0.76%  |   0.03%  |  13.32%  |  72.93%
--------------------------------------------------