Coath, Christopher D., Robert CJ Steele, and W. Fred Lunnon. "Statistical bias in isotope ratios." Journal of Analytical Atomic Spectrometry 28.1 (2013): 52-58.
$\frac{X}{Y+1}$ has the expectation value of $(\lambda_1/\lambda_2)(1-\exp{-\lambda_2})$ NOT CONVINCED THIS IS TRUE
In [68]:
import pymc3 as pm
import numpy as np
import seaborn as nsn
import pandas as pd
import matplotlib.pyplot as plt
import scipy
import seaborn as sns
In [69]:
rate1 = 15
rate2 = rate1/11
print(rate1, rate2, rate1/rate2, (rate1/rate2)*(1-np.exp(-rate2)))
d1 = np.random.poisson(rate1, size=1000)
d2 = np.random.poisson(rate2, size=1000)
In [80]:
with pm.Model() as model1:
r = pm.Uniform('rate', 0, 100, shape=2)
p = pm.Poisson('data', r, shape=2, observed=np.vstack((d1, d2)).T)
ratio = pm.Deterministic('ratio', r[0]/(r[1]+1))
sum1 = pm.Deterministic('sum1', pm.math.sum(r[0]))
sum2 = pm.Deterministic('sum2', pm.math.sum(r[1]))
start = pm.find_MAP()
trace1 = pm.sample(10000, start=start, njobs=5)
In [81]:
pm.traceplot(trace1, combined=True)
Out[81]:
In [82]:
pm.summary(trace1)
Out[82]:
In [83]:
(trace1['rate'][:,0].mean(),
trace1['rate'][:,1].mean(),
(trace1['rate'][:,0].mean()/trace1['rate'][:,1].mean()),
(1-np.exp(-trace1['rate'][:,1].mean())),
(trace1['rate'][:,0].mean()/trace1['rate'][:,1].mean()) * (1-np.exp(-trace1['rate'][:,1].mean())) )
Out[83]:
In [84]:
sns.distplot(trace1['ratio'])
plt.axvline( (trace1['rate'][:,0].mean()/trace1['rate'][:,1].mean()) * (1-np.exp(-trace1['rate'][:,1].mean())) )
Out[84]:
In [85]:
sns.distplot(trace1['sum1']/trace1['sum2'])
plt.axvline(11)
Out[85]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [6]:
rate1 = 15*np.sin(np.linspace(0, np.pi, 1000))+4
rate2 = rate1/11
d1 = scipy.random.poisson(rate1)
d2 = scipy.random.poisson(rate2)
plt.plot(d1)
plt.plot(d2)
Out[6]:
In [7]:
plt.plot(d1/d2)
Out[7]:
In [8]:
plt.hist((d1/d2)[np.isfinite(d1/d2)], 20)
Out[8]:
In [87]:
with pm.Model() as model2:
r1 = pm.Uniform('rate1', 0, 100, shape=len(d1))
r2 = pm.Uniform('rate2', 0, 100, shape=len(d2))
p1 = pm.Poisson('data1', r1, shape=len(d1), observed=d1)
p2 = pm.Poisson('data2', r2, shape=len(d2), observed=d2)
ratio = pm.Deterministic('ratio', r1/(r2+1))
sum1 = pm.Deterministic('sum1', pm.math.sum(r1))
sum2 = pm.Deterministic('sum2', pm.math.sum(r2))
start = pm.find_MAP()
trace2 = pm.sample(10000, start=start, njobs=5)
In [88]:
pm.summary(trace2, varnames=['ratio'])
Out[88]:
In [90]:
sns.distplot(trace2['sum1'].flatten()/trace2['sum2'].flatten())
plt.axvline(11)
Out[90]:
In [ ]: