In [1]:
import datetime
print(datetime.datetime.now().isoformat())
Using a statistical model of the relationship between the slot filling at 464keV vs Dst.
We model the number of slot filling events as a random sample from a binomial distribution meaning that there is a probability of the event occuring and there may be a parameter that changes the probability. This is not meant as a correlation but as a success failure creiteria.
From Wikipedia: In probability theory and statistics, the binomial distribution with parameters n and p is the discrete probability distribution of the number of successes in a sequence of n independent yes/no experiments, each of which yields success with probability p.
$ change \sim Bin(n,p) \\ logit(p) = \alpha + \beta x \\ a \sim N(0,5) \\ \beta \sim N(0,10) $
where we set vague priors for $\alpha$ and $\beta$, the parameters for the logistic model.
This is the same technique used in the estimation of deaths due to a concentration of a chemical.
In [2]:
# http://onlinelibrary.wiley.com/doi/10.1002/2016JA022652/epdf
import pymc
from pprint import pprint
import numpy as np
import matplotlib.pyplot as plt
import spacepy.plot as spp
Data delivered by Geoff Reeves 9/6/2016
In [3]:
# min_Dst, min_L
data = np.asarray([
# data from 20160906 email from Geoff
# for these slot filling is L=2.3, set to 0 below
# 65.000, 3.8000,
# 50.000, 3.7000,
# 67.000, 3.5000,
# 61.000, 3.4000,
# 77.000, 3.2000,
# 99.000, 2.8900,
# 87.000, 2.8000,
# 98.000, 2.8000,
# 96.000, 2.8000,
# 93.000, 2.3000,
# 92.000, 2.3000,
# 225.00, 2.3000,
# 206.00, 2.3000,
# 125.00, 2.3000,
# data from 20160907 email form Geoff
# for these slot filling is L=0
53.000, 4.0000,
54.000, 3.8000,
133.00, 0.0000,
61.000, 3.6000,
66.000, 3.3000,
56.000, 3.7000,
118.00, 2.9000,
73.000, 3.2000,
96.000, 2.8000,
72.000, 3.5000,
44.000, 4.0000,
54.000, 3.8000,
67.000, 3.1000,
80.000, 3.7000,
66.000, 4.3000,
116.00, 3.2000,
94.000, 3.4000,
81.000, 3.6000,
65.000, 4.0000,
38.000, 3.7000,
80.000, 3.7000,
76.000, 3.5000,
71.000, 3.5000,
99.000, 2.9000,
47.000, 3.7000,
54.000, 3.3000,
224.00, 0.0000,
76.000, 3.2000,
75.000, 3.2000,
73.000, 3.6000,
203.00, 0.0000,
68.000, 3.5000,
61.000, 3.4000,
63.000, 3.4000,
84.000, 2.9000,
92.000, 0.0000,
97.000, 0.0000,
75.000, 3.2000,
123.00, 0.0000,
97.000, 2.7000,
60.000, 3.7000,
171.00, 0.0000,
111.00, 3.3000,
93.000, 3.1000,
57.000, 3.5000,
97.000, 3.1000,
55.000, 3.8000,
30.000, 3.8000,
83.000, 0.0000,
]).reshape((-1,2))
dst = data[:,0]
minL = data[:,1]
print(dst, minL, data.dtype)
# ind_fill = minL <= 2.3
# minL[ind_fill] = 0 # set them to zero
In [4]:
# make bins in Dst
dst_bins = np.arange(25, 300, 1)
print(dst_bins)
dst_bins_centers = np.asarray([dst_bins[:-1] + np.diff(dst_bins)/2]).T[:,0]
print(dst_bins_centers, dst_bins_centers.shape)
n_events_dig = np.digitize(dst, dst_bins)
print(n_events_dig)
n_events = np.zeros(len(dst_bins)-1)
success = np.zeros_like(n_events)
for i, v in enumerate(np.unique(n_events_dig)):
n_events[v-1] = np.sum(n_events_dig==v)
success[v-1] = np.sum(minL[n_events_dig==v] <= 2.3)
print(n_events)
print(success)
In [5]:
plt.plot(dst, minL, 'o')
plt.xlim((25, 250))
plt.ylim((-0.5, 5.0))
plt.xlabel('Min Dst')
plt.ylabel('Min L Shell')
Out[5]:
In [6]:
plt.plot(dst, minL, 'o')
plt.xlim((25, 250))
plt.ylim((-0.5, 5.0))
plt.xlabel('Min Dst')
plt.ylabel('Min L Shell')
# for v in dst_bins:
# plt.axvline(v, lw=0.5)
plt.axhline(2.4, c='r', lw=1)
Out[6]:
In [7]:
# define invlogit function
def invlogit(x):
return pymc.exp(x) / (1 + pymc.exp(x))
Setup the Bayesian model accorind to the description above. Run the Markov chain monte carlo (MCMC) to sample the posterior distributions for $\alpha$ and $\beta$
In [8]:
# define priors
# these are wide uninformative priors
# alpha = pymc.Normal('alpha', mu=0, tau=1.0/5**2)
# beta = pymc.Normal('beta', mu=0, tau=1.0/10**2)
alpha = pymc.Uniform('alpha', -500, 100, value = 1e-2)
beta = pymc.Uniform('beta', -100, 100, value = 1e-2)
# cannot feed in zero events
ind = n_events > 0
print(n_events.shape, dst_bins_centers.shape, success.shape, )
# define likelihood
p = pymc.InvLogit('p', alpha + beta*dst_bins_centers[ind])
print('n_events', n_events[ind], 'p', p.value, 'success', success[ind], )
y = pymc.Binomial('y_obs', n=n_events[ind], p=p, value=success[ind], observed=True)
In [9]:
# inference
m = pymc.Model([alpha, beta, y])
mc = pymc.MCMC(m)
mc.sample(iter=500000, burn=1000, burn_till_tuned=True, thin=110)
Make diagnostic plots of the posteriour distributions as created using MCMC.
In [10]:
pymc.Matplot.plot(mc)
In [11]:
pprint(mc.stats())
In [12]:
xp = np.linspace(dst_bins_centers[ind].min(), dst_bins_centers[ind].max(), 100)
a = alpha.stats()['mean']
b = beta.stats()['mean']
y_val = invlogit(a + b*xp).value
plt.plot(xp, y_val)
plt.scatter(dst_bins_centers[ind], success[ind]/n_events[ind], s=50);
plt.xlabel('Minimum Dst')
plt.ylabel('Probability slot is filled')
plt.gca().ticklabel_format(useOffset=False)
In [13]:
# get the minimum Dst where 99% should be successes
for percentage in [50,75,90,95,99]:
ind99 = y_val >= percentage/100
minDst99 = xp[ind99][0]
print('At a minimum Dst of {0:0.0f}nT it is predicted to have a {1}% percent of a slot filling'.format(minDst99, percentage))
The 20160906 data gave these results:
At a minimum Dst of 100nT it is predicted to have a 50% percent of a slot filling
At a minimum Dst of 105nT it is predicted to have a 75% percent of a slot filling
At a minimum Dst of 111nT it is predicted to have a 90% percent of a slot filling
At a minimum Dst of 114nT it is predicted to have a 95% percent of a slot filling
At a minimum Dst of 123nT it is predicted to have a 99% percent of a slot filling
In [ ]:
In [14]:
# one should be able to get estimates of the line uncertainity
ilu = np.empty((1000, len(xp)), dtype=float)
for ii, v in enumerate(np.random.random_integers(0, len(alpha.trace[:])-1, 1000)):
ilu[ii] = invlogit(alpha.trace[:][v] + beta.trace[:][v]*xp).value
In [15]:
xp = np.linspace(dst_bins_centers[ind].min(), dst_bins_centers[ind].max(), 100)
for v in ilu:
plt.plot(xp, v, alpha=.01, c='r')
a = alpha.stats()['mean']
b = beta.stats()['mean']
plt.plot(xp, invlogit(a + b*xp).value)
a = alpha.stats()['quantiles'][50]
b = beta.stats()['quantiles'][50]
plt.plot(xp, invlogit(a + b*xp).value, c='y')
plt.scatter(dst_bins_centers[ind], success[ind]/n_events[ind], s=50);
plt.xlabel('Minimum Dst')
plt.ylabel('Probability slot is filled')
Out[15]:
Same as previous figure with the red lines overlayed as 100 joint draws from the model posterior in order to show spread.
In [16]:
# make bins in Dst
dst_bins = np.arange(25, 300, 10)
print(dst_bins)
dst_bins_centers = np.asarray([dst_bins[:-1] + np.diff(dst_bins)/2]).T[:,0]
print(dst_bins_centers, dst_bins_centers.shape)
n_events_dig = np.digitize(dst, dst_bins)
print(n_events_dig)
n_events = np.zeros(len(dst_bins)-1)
success = np.zeros_like(n_events)
for i, v in enumerate(np.unique(n_events_dig)):
n_events[v-1] = np.sum(n_events_dig==v)
success[v-1] = np.sum(minL[n_events_dig==v] <= 2.3)
print(n_events)
print(success)
In [17]:
plt.plot(dst, minL, 'o')
plt.xlim((25, 250))
plt.ylim((-0.5, 5.0))
plt.xlabel('Min Dst')
plt.ylabel('Min L Shell')
for v in dst_bins:
plt.axvline(v, lw=0.5)
plt.axhline(2.4, c='r', lw=1)
Out[17]:
In [18]:
# define priors
# these are wide uninformative priors
# alpha = pymc.Normal('alpha', mu=0, tau=1.0/5**2)
# beta = pymc.Normal('beta', mu=0, tau=1.0/10**2)
alpha = pymc.Uniform('alpha', -500, 100, value = 1e-2)
beta = pymc.Uniform('beta', -100, 100, value = 1e-2)
# cannot feed in zero events
ind = n_events > 0
print(n_events.shape, dst_bins_centers.shape, success.shape, )
# define likelihood
p = pymc.InvLogit('p', alpha + beta*dst_bins_centers[ind])
print('n_events', n_events[ind], 'p', p.value, 'success', success[ind], )
y = pymc.Binomial('y_obs', n=n_events[ind], p=p, value=success[ind], observed=True)
In [19]:
# inference
m = pymc.Model([alpha, beta, y])
mc = pymc.MCMC(m)
mc.sample(iter=500000, burn=1000, burn_till_tuned=True, thin=110)
In [20]:
pymc.Matplot.plot(mc)
pprint(mc.stats())
In [21]:
xp = np.linspace(dst_bins_centers[ind].min(), dst_bins_centers[ind].max(), 100)
a = alpha.stats()['mean']
b = beta.stats()['mean']
y_val = invlogit(a + b*xp).value
plt.plot(xp, y_val)
plt.scatter(dst_bins_centers[ind], success[ind]/n_events[ind], s=50);
plt.xlabel('Minimum Dst')
plt.ylabel('Probability slot is filled')
plt.gca().ticklabel_format(useOffset=False)
In [22]:
# get the minimum Dst where 99% should be successes
for percentage in [1, 25, 50,75,90,95,99]:
ind99 = y_val >= percentage/100
minDst99 = xp[ind99][0]
print('At a minimum Dst of {0:0.0f}nT it is predicted to have a {1}% percent of a slot filling'.format(minDst99, percentage))
In [23]:
# one should be able to get estimates of the line uncertainity
ilu = np.empty((1000, len(xp)), dtype=float)
for ii, v in enumerate(np.random.random_integers(0, len(alpha.trace[:])-1, 1000)):
ilu[ii] = invlogit(alpha.trace[:][v] + beta.trace[:][v]*xp).value
In [24]:
xp = np.linspace(dst_bins_centers[ind].min(), dst_bins_centers[ind].max(), 100)
for v in ilu:
plt.plot(xp, v, alpha=.01, c='r')
a = alpha.stats()['mean']
b = beta.stats()['mean']
plt.plot(xp, invlogit(a + b*xp).value)
a = alpha.stats()['quantiles'][50]
b = beta.stats()['quantiles'][50]
plt.plot(xp, invlogit(a + b*xp).value, c='y')
plt.scatter(dst_bins_centers[ind], success[ind]/n_events[ind], s=50);
plt.xlabel('Minimum Dst')
plt.ylabel('Probability slot is filled')
Out[24]:
In [ ]:
In [25]:
# make bins in Dst
dst_bins = np.arange(25, 300, 25)
print(dst_bins)
dst_bins_centers = np.asarray([dst_bins[:-1] + np.diff(dst_bins)/2]).T[:,0]
print(dst_bins_centers, dst_bins_centers.shape)
n_events_dig = np.digitize(dst, dst_bins)
print(n_events_dig)
n_events = np.zeros(len(dst_bins)-1)
success = np.zeros_like(n_events)
for i, v in enumerate(np.unique(n_events_dig)):
n_events[v-1] = np.sum(n_events_dig==v)
success[v-1] = np.sum(minL[n_events_dig==v] <= 2.3)
print(n_events)
print(success)
In [26]:
plt.plot(dst, minL, 'o')
plt.xlim((25, 250))
plt.ylim((-0.5, 5.0))
plt.xlabel('Min Dst')
plt.ylabel('Min L Shell')
for v in dst_bins:
plt.axvline(v, lw=0.5)
plt.axhline(2.4, c='r', lw=1)
Out[26]:
In [27]:
# define priors
# these are wide uninformative priors
# alpha = pymc.Normal('alpha', mu=0, tau=1.0/5**2)
# beta = pymc.Normal('beta', mu=0, tau=1.0/10**2)
alpha = pymc.Uniform('alpha', -500, 100, value = 1e-2)
beta = pymc.Uniform('beta', -100, 100, value = 1e-2)
# cannot feed in zero events
ind = n_events > 0
print(n_events.shape, dst_bins_centers.shape, success.shape, )
# define likelihood
p = pymc.InvLogit('p', alpha + beta*dst_bins_centers[ind])
print('n_events', n_events[ind], 'p', p.value, 'success', success[ind], )
y = pymc.Binomial('y_obs', n=n_events[ind], p=p, value=success[ind], observed=True)
In [28]:
# inference
m = pymc.Model([alpha, beta, y])
mc = pymc.MCMC(m)
mc.sample(iter=500000, burn=1000, burn_till_tuned=True, thin=110)
In [29]:
pymc.Matplot.plot(mc)
pprint(mc.stats())
In [30]:
xp = np.linspace(dst_bins_centers[ind].min(), dst_bins_centers[ind].max(), 100)
a = alpha.stats()['mean']
b = beta.stats()['mean']
y_val = invlogit(a + b*xp).value
plt.plot(xp, y_val)
plt.scatter(dst_bins_centers[ind], success[ind]/n_events[ind], s=50);
plt.xlabel('Minimum Dst')
plt.ylabel('Probability slot is filled')
plt.gca().ticklabel_format(useOffset=False)
In [31]:
# get the minimum Dst where 99% should be successes
for percentage in [1, 25, 50,75,90,95,99]:
ind99 = y_val >= percentage/100
minDst99 = xp[ind99][0]
print('At a minimum Dst of {0:0.0f}nT it is predicted to have a {1}% percent of a slot filling'.format(minDst99, percentage))
In [32]:
# one should be able to get estimates of the line uncertainity
ilu = np.empty((1000, len(xp)), dtype=float)
for ii, v in enumerate(np.random.random_integers(0, len(alpha.trace[:])-1, 1000)):
ilu[ii] = invlogit(alpha.trace[:][v] + beta.trace[:][v]*xp).value
In [33]:
xp = np.linspace(dst_bins_centers[ind].min(), dst_bins_centers[ind].max(), 100)
for v in ilu:
plt.plot(xp, v, alpha=.01, c='r')
a = alpha.stats()['mean']
b = beta.stats()['mean']
plt.plot(xp, invlogit(a + b*xp).value)
a = alpha.stats()['quantiles'][50]
b = beta.stats()['quantiles'][50]
plt.plot(xp, invlogit(a + b*xp).value, c='y')
plt.scatter(dst_bins_centers[ind], success[ind]/n_events[ind], s=50);
plt.xlabel('Minimum Dst')
plt.ylabel('Probability slot is filled')
Out[33]: