Bi 1x: Maximum likelihood estimation¶

This tutorial was generated from a Jupyter notebook, which can be downloaded here

In this tutorial, we will use our list of bead loss times to get an estimate for the characteristic cleavage time by RAG (for the concentration of our given experiment). We begin by importing our usual modules.

In [1]:

import numpy as np
import scipy.optimize

import iqplot

import bokeh.io
bokeh.io.output_notebook()

BokehJS 3.1.1 successfully loaded.

And, of course, we need to load in the data set we saved from the Image Processing III tutorial. In the event that you didn't complete that tutorial, you can download the the bead loss times here.

In [2]:

# Data from image processing
t_lost = np.loadtxt('t_lost_2016.csv')
m = int(np.loadtxt('number_of_uncut_2016.csv'))
t = np.loadtxt('time_points_2016.csv')
n = len(t_lost)
t_end = t[-1]

The distribution of observed cleavage times¶

We assume that cleavage is a Poisson process. A Poisson process is a sequence of events (called arrivals) that occur over time such that the amount of time between any two arrivals is independent of all other inter-arrival times. This means that cleavage is a stochastic event that has no memory, nor knowledge of the future.

This is essentially our physical model for cleavage. We do not know the molecular details, only that the cleavage events are rare and random.

The amount of time to wait for an event in a Poisson process is exponentially distributed. We will not derive that fact here, but just state it. That means that the probability of observing a bead being lost between time $t$ and $t + \mathrm{d}t$ , where $\mathrm{d}t$ is a differential time difference, is $\mathrm{d}t\,f(t;\beta)$ , where

\begin{aligned} (1) & f (t; β) = β e^{- β t} . \end{aligned}

$\begin{align} f(t ; \beta) = \beta\,\mathrm{e}^{-\beta t}. \end{align}$

The function $f(t;\beta)$ is called a probability density function, or PDF, and has units of inverse time. The parameter $\beta$ , which has units of inverse time, and characterizes how rapid cleavage is. It is proportional to the chemical rate constant from mass action kinetics that you may be accustomed to. Specifically, the constant of proportionality contains the concentrations of enzyme, DNA, and volume of the system.

The problem with this argument is that beads may be lost to processes other than cleavage by RAG protein. We do a control experiment without RAG to estimate the rate of bead loss in the absence of RAG, the rate constant of which we will define as $\beta_\mathrm{auto}$ . So, now we have two Poisson processes that could arrive and a bead could be lost. I.e., the probability density function for bead loss due to non-RAG-mediated bead loss is

\begin{aligned} (2) & f (t; β_{a u t o}) = β_{a u t o} e^{- β_{a u t o} t} . \end{aligned}

$\begin{align} f(t ; \beta_\mathrm{auto}) = \beta_\mathrm{auto}\,\mathrm{e}^{-\beta_\mathrm{auto}t}. \end{align}$

Since we are observing loss, possibly due to either process, we really want is the probability distribution of bead loss due to either RAG cleavage or other processes. This is a more general problem: how is the waiting time for the arrival of one of two different Poisson processes distributed? To derive this, we first write down the probability that a Poisson process has not arrived before time $t'$ . That is, we want $P(t > t'; \beta)$ .

\begin{aligned} (3) & P (t > t^{'}; β) = \int_{t^{'}}^{\infty} d t f (t; β) = \int_{t^{'}}^{\infty} d t β e^{- β t} = e^{- β t^{'}} . \end{aligned}

$\begin{align} P(t > t'; \beta) = \int_{t'}^\infty \mathrm{d}t\,f(t;\beta) = \int_{t'}^\infty \mathrm{d}t\,\beta\mathrm{e}^{-\beta t} = \mathrm{e}^{-\beta t'}. \end{align}$

Now, if the two processes, which we will label with $1$ and $2$ , are independent of each other, the probability that neither have arrived before time $t'$ is

\begin{aligned} (4) & P (t_{1}, t_{2} > t^{'}; β_{1}, β_{2}) & = P (t_{1} > t^{'}; β_{1}) P (t_{2} > t^{'}; β_{2}) = e^{- β_{1} t^{'}} e^{- β_{2} t^{'}} = e^{- (β_{1} + β_{2}) t^{'}} . \end{aligned}

$\begin{align} P(t_1, t_2 > t'; \beta_1, \beta_2) &= P(t_1>t'; \beta_1)\,P(t_2>t'; \beta_2) = \mathrm{e}^{-\beta_1 t'}\,\mathrm{e}^{-\beta_2 t'} = \mathrm{e}^{-\left(\beta_1 + \beta_2\right)t'}. \end{align}$

But this is the expression for the probability that a single Poisson process with rate $\beta = \beta_1 + \beta_2$ does not arrive by time $t'$ . So, the distribution for the arrival time of either Poisson process in time $t$ is that of a single Poisson process with time constant $\beta = \beta_1 + \beta_2$ . So, we have

\begin{aligned} (5) & f (t; β, β_{a u t o}) = (β + β_{a u t o}) e^{- (β_{1} + β_{2}) t} . \end{aligned}

$\begin{align} f(t ; \beta, \beta_\mathrm{auto}) = \left(\beta + \beta_\mathrm{auto}\right)\mathrm{e}^{-\left(\beta_1 + \beta_2\right)t}. \end{align}$

For convenience going forward, we will define $\beta$ to contain both the rate of cleavage and that of auto-detachment, knowing that we can subtract off the rate of auto-detachment obtained from the control experiment at the end of the analysis. So, we will use

\begin{aligned} (6) & f (t; β) = β e^{- β t} . \end{aligned}

$\begin{align} f(t; \beta) = \beta \mathrm{e}^{-\beta t}. \end{align}$

If we observe an exponential distribution of cleavage times, we have strong evidence that cleavage is well-described by a Poisson process; it is just a stochastic memoryless process. If, however, the cleavage times are not exponentially distributed, we have to think hard about what might be at play in the dynamics of RAG-mediated cleavage.

We have the PDF describing bead loss at time $t$ , but we want the probability of our observed bead loss time. We need to be more precise in what we mean here. We take images at time points $t_0$ , $t_1$ , ..., $t_\mathrm{end}$ . We can observe a bead is missing at any of these times points except $t_0$ . Let $t_i$ be the time where we observe a given bead is missing. What is the probability of this observation? It is the probability that the bead was lost between time $t_{i-1}$ and time $t_i$ .

\begin{aligned} (7) & P (observed loss at t_{i}; β) = \int_{t_{i - 1}}^{t_{i}} d t^{'} f (t^{'}; β) = \int_{t_{i - 1}}^{t_{i}} d t^{'} β e^{- β t^{'}} = e^{- β t_{i - 1}} - e^{- β t_{i}} = e^{- β t_{i}} (e^{β (t_{i} - t_{i - 1})} - 1) . \end{aligned}

$\begin{align} P(\text{observed loss at } t_i ; \beta) = \int_{t_{i-1}}^{t_i}\mathrm{d}t'\,f(t';\beta) = \int_{t_{i-1}}^{t_i}\mathrm{d}t'\,\beta \mathrm{e}^{-\beta t'} = \mathrm{e}^{-\beta t_{i-1}} - \mathrm{e}^{-\beta t_i} = \mathrm{e}^{-\beta t_i}\left(\mathrm{e}^{\beta (t_i - t_{i-1})} - 1\right). \end{align}$

If we have the same time interval $\Delta t$ between each image we take, then $t_i - t_{i-1} = \Delta t$ for all $i$ , and this expression is

\begin{aligned} (8) & P (observed loss at t_{i}; β) = e^{- β t_{i}} (e^{β Δ t} - 1) . \end{aligned}

$\begin{align} P(\text{observed loss at } t_i ; \beta) = \mathrm{e}^{-\beta t_i}\left(\mathrm{e}^{\beta \Delta t} - 1\right). \end{align}$

Note that we may not always have the same time interval between frames because of loss of focus or other issues, we will make this assumption going forward.

Since the bead loss events are all independent, we can write down the probability of observing all of our lost beads as the product of these probabilities. Let $t_{ij}$ be the time at which we observe bead $j$ going missing. Let $\mathbf{t}$ be a list of the time points at which we observe each of of the $n$ beads we lost. Then,

\begin{aligned} (9) & P (t; β) = \prod_{j = 1}^{n} e^{- β t_{i j}} (e^{β Δ t_{i j}} - 1) = {(e^{β Δ t} - 1)}^{n} \prod_{j = 1}^{n} e^{- β t_{i j}}, \end{aligned}

$\begin{align} P(\mathbf{t} ; \beta) = \prod_{j=1}^n\mathrm{e}^{-\beta t_{ij}}\left(\mathrm{e}^{\beta \Delta t_{ij}} - 1\right) = \left(\mathrm{e}^{\beta \Delta t} - 1\right)^n\prod_{j=1}^n\mathrm{e}^{-\beta t_{ij}}, \end{align}$

where $\Delta t_{ij} = t_{ij} - t_{i-1,j}$ is the time interval between the image where bead $j$ is first absent and the last image where it was present. The last equality in the above expression holds if we have constant time intervals between images. We note that

\begin{aligned} (10) & \sum_{j = 1}^{n} t_{i j} = n {\bar{t}}_{i}, \end{aligned}

$\begin{align} \sum_{j=1}^n t_{ij} = n \bar{t}_i, \end{align}$

where $\bar{t}_i$ is the average of the times where we observed bead loss. As a result,

\begin{aligned} (11) & P (t; β) = e^{- n β {\bar{t}}_{i}} \prod_{j = 1}^{n} (e^{β Δ t_{i j}} - 1) = e^{- n β {\bar{t}}_{i}} {(e^{β Δ t} - 1)}^{n}, \end{aligned}

$\begin{align} P(\mathbf{t} ; \beta) = \mathrm{e}^{-n \beta \bar{t}_i}\prod_{j=1}^n\left(\mathrm{e}^{\beta \Delta t_{ij}} - 1\right) = \mathrm{e}^{-n \beta \bar{t}_i}\left(\mathrm{e}^{\beta \Delta t} - 1\right)^n, \end{align}$

where again the last equality holds if we have constant time intervals between images.

Some beads are not cleaved¶

We only run our movie for time $t_\mathrm{end}$ . (In your experiment, we typically set $t_\mathrm{end} = 3$ hours.) So there may be beads that could be cleaved if we waited a while. These are data, too! The probability that a given bead is not cleaved before $t_\mathrm{end}$ , which means that we observe it intact in the final frame of the movie, is

\begin{aligned} (12) & P (t > t_{e n d}; β) = \int_{t_{e n d}}^{\infty} d t^{'} f (t^{'}; β) = \int_{t_{e n d}}^{\infty} d t^{'} β e^{- β t^{'}} = e^{- β t_{e n d}} . \end{aligned}

$\begin{align} P(t > t_\mathrm{end} ; \beta) = \int_{t_\mathrm{end}}^\infty \mathrm{d}t'\,f(t';\beta) = \int_{t_\mathrm{end}}^\infty \mathrm{d}t'\,\beta \mathrm{e}^{-\beta t'} = \mathrm{e}^{-\beta t_\mathrm{end}}. \end{align}$

If there are $m$ such beads, then the probability of observing our $n$ cleaved beads and our $m$ uncleaved beads in the experiment is

\begin{aligned} (13) & P (t, m; β) = {(P (t > t_{e n d}; β))}^{m} P (t; β) = e^{- m β t_{e n d}} e^{- n β {\bar{t}}_{i}} \prod_{j = 1}^{n} (e^{β Δ t_{i j}} - 1) = e^{- m β t_{e n d}} e^{- n β {\bar{t}}_{i}} {(e^{β Δ t} - 1)}^{n}, \end{aligned}

$\begin{align} P(\mathbf{t},m;\beta) = \left(P(t > t_\mathrm{end} ; \beta)\right)^mP(\mathbf{t} ; \beta) = \mathrm{e}^{-m\beta t_\mathrm{end}}\mathrm{e}^{-n \beta \bar{t}_i}\prod_{j=1}^n\left(\mathrm{e}^{\beta \Delta t_{ij}} - 1\right) = \mathrm{e}^{-m\beta t_\mathrm{end}}\mathrm{e}^{-n \beta \bar{t}_i}\left(\mathrm{e}^{\beta \Delta t} - 1\right)^n, \end{align}$

where again the last equality holds if we have constant time intervals between images.

Maximum likelihood estimation¶

To get an estimate for $\beta$ , the rate of cleavage (which, recall, is actually the rate of cleavage plus the rate of non-enzyme-mediated bead loss), we choose the value of $\beta$ such that the expression for $P(\mathbf{t},m;\beta)$ is maximized. This is called a maximum likelihood estimate.

The value of $\beta$ that maximizes $P(\mathbf{t},m;\beta)$ is the same that maximized $\ln P(\mathbf{t},m;\beta)$ , and the logarithm is much easier to work with.

\begin{aligned} (14) & \ln P (t, m; β) = - m β t_{e n d} - n β {\bar{t}}_{i} + \sum_{j = 1}^{n} \ln (e^{β Δ t_{i j}} - 1) . \end{aligned}

$\begin{align} \ln P(\mathbf{t},m;\beta) = -m\beta t_\mathrm{end} - n \beta \bar{t}_i + \sum_{j=1}^n \ln \left(\mathrm{e}^{\beta \Delta t_{ij}} - 1\right). \end{align}$

This is maximal where its derivative vanishes, that is where

\begin{aligned} (15) & \frac{d}{d β} \ln P (t, m; β) = - m t_{e n d} - n {\bar{t}}_{i} + \sum_{j = 1}^{n} \frac{Δ t_{i j} e^{β Δ t_{i j}}}{e^{β Δ t_{i j}} - 1} = 0. \end{aligned}

$\begin{align} \frac{\mathrm{d}}{\mathrm{d}\beta}\,\ln P(\mathbf{t},m;\beta) = -m t_\mathrm{end} - n \bar{t}_i + \sum_{j=1}^n \frac{\Delta t_{ij} \mathrm{e}^{\beta \Delta t_{ij}}}{\mathrm{e}^{\beta \Delta t_{ij}} - 1} = 0. \end{align}$

In the case where we have evenly spaced time points, this is

\begin{aligned} (16) & \frac{d}{d β} \ln P (t, m; β) = - m t_{e n d} - n {\bar{t}}_{i} + \frac{n Δ t e^{β Δ t}}{e^{β Δ t} - 1} = 0. \end{aligned}

$\begin{align} \frac{\mathrm{d}}{\mathrm{d}\beta}\,\ln P(\mathbf{t},m;\beta) = -m t_\mathrm{end} - n \bar{t}_i + \frac{n\Delta t\mathrm{e}^{\beta \Delta t}}{\mathrm{e}^{\beta \Delta t} - 1} = 0. \end{align}$

Approximate MLE¶

To solve for the MLE, we need find the value of $\beta$ for which

\begin{aligned} (17) & - m t_{e n d} - n {\bar{t}}_{i} + \sum_{j = 1}^{n} \frac{Δ t_{i j} e^{β Δ t_{i j}}}{e^{β Δ t_{i j}} - 1} = 0. \end{aligned}

$\begin{align} -m t_\mathrm{end} - n \bar{t}_i + \sum_{j=1}^n \frac{\Delta t_{ij} \mathrm{e}^{\beta \Delta t_{ij}}}{\mathrm{e}^{\beta \Delta t_{ij}} - 1} = 0. \end{align}$

We now consider an approximate solution for the special case where we have evenly spaced time points so that $\Delta t_{ij} = \Delta t$ and rapid image acquisition so that $\beta \Delta t \ll 1$ . In this case,

\begin{aligned} (18) & e^{β Δ t} \approx 1 + β Δ t, \end{aligned}

$\begin{align} \mathrm{e}^{\beta \Delta t} \approx 1 + \beta \Delta t, \end{align}$

so that

\begin{aligned} (19) & - m t_{e n d} - n {\bar{t}}_{i} + \frac{n Δ t (1 + β Δ t)}{1 + β Δ t - 1} \approx - m t_{e n d} - n {\bar{t}}_{i} + \frac{n}{β} \approx 0. \end{aligned}

$\begin{align} -m t_\mathrm{end} - n \bar{t}_i + \frac{n\Delta t(1 + \beta \Delta t)}{1 + \beta \Delta t - 1} \approx -m t_\mathrm{end} - n \bar{t}_i + \frac{n}{\beta} \approx 0. \end{align}$

Solving for $\beta$ gives

\begin{aligned} (20) & β = {({\bar{t}}_{i} + \frac{m}{n} t_{e n d})}^{- 1} . \end{aligned}

$\begin{align} \beta = \left(\bar{t}_i + \frac{m}{n}\,t_\mathrm{end}\right)^{-1}. \end{align}$

So, if we waited a long time such that we had essentially no beads remaining, the $\beta \approx 1/\bar{t}_i$ . In other words, the rate $\beta$ is best estimated by the inverse of the average time of cleavage.

We can compute the approximate MLE for $\beta$ using the data.

In [3]:

beta_mle_approx = 1 / (np.mean(t_lost) + m / n * t_end)

print(beta_mle_approx)

0.013887675779544067

The characteristic cutting time is $1 / \beta_\mathrm{MLE}$ , which we now compute.

In [4]:

1 / beta_mle_approx

Out[4]:

72.0062893081761

The characteristic cut time is about 72 minutes.

Note that in this approximation, we effectively assumed that even though we have discrete sampling times, we are able to measure the exact time of losing a bead. This is clear by forgoing integrating over the interval from time $t_{i-1}$ to $t_i$ for a bead loss observed at time $t_i$ , and instead directly using the probability density function such that

\begin{aligned} (21) & P (t, m; β) = e^{- m β t_{e n d}} \prod_{i = j}^{n} β e^{- β t_{j}} = β^{n} e^{- m β t_{e n d}} e^{- n β \bar{t}} . \end{aligned}

$\begin{align} P(\mathbf{t},m;\beta) = \mathrm{e}^{-m\beta t_\mathrm{end}} \prod_{i=j}^n \beta \mathrm{e}^{-\beta t_j} = \beta^n\,\mathrm{e}^{-m\beta t_\mathrm{end}}\,\mathrm{e}^{-n\beta \bar{t}}. \end{align}$

The logarithm of this is

\begin{aligned} (22) & \ln P (t, m; β) = - m β t_{e n d} - n β \bar{t} + n \ln β . \end{aligned}

$\begin{align} \ln P(\mathbf{t},m;\beta) = -m\beta t_\mathrm{end} - n\beta \bar{t} + n\ln \beta. \end{align}$

Differentiating with respect to $\beta$ and setting to zero gives

\begin{aligned} (23) & \frac{d}{d β} = - m t_{e n d} - n \bar{t} + \frac{n}{β} = 0. \end{aligned}

$\begin{align} \frac{\mathrm{d}}{\mathrm{d}\beta} = -mt_\mathrm{end} - n\bar{t} + \frac{n}{\beta} = 0. \end{align}$

Solving for $\beta$ gives

\begin{aligned} (24) & β = {({\bar{t}}_{i} + \frac{m}{n} t_{e n d})}^{- 1}, \end{aligned}

$\begin{align} \beta = \left(\bar{t}_i + \frac{m}{n}\,t_\mathrm{end}\right)^{-1}, \end{align}$

the same result as formally making the approximation above.

Computing the MLE¶

We now demonstrate how to compute the MLE exactly. Our task is to write a function that computes the function of $\beta$ we are trying to find the root of, namely

\begin{aligned} (25) & g (β) = - m t_{e n d} - n {\bar{t}}_{i} + \sum_{j = 1}^{n} \frac{Δ t_{i j} e^{β Δ t_{i j}}}{e^{β Δ t_{i j}} - 1} . \end{aligned}

$\begin{align} g(\beta) = -m t_\mathrm{end} - n \bar{t}_i + \sum_{j=1}^n \frac{\Delta t_{ij} \mathrm{e}^{\beta \Delta t_{ij}}}{\mathrm{e}^{\beta \Delta t_{ij}} - 1}. \end{align}$

We first compute some convenient quantities from the data.

In [5]:

# Array of Δt's
delta_t = np.concatenate(((0,), np.diff(t)))

# Index of time points where cut
i = np.array([np.where(t == tl)[0][0] for tl in t_lost])

# The Δtij's
delta_t_ij = delta_t[i]

# Mean observation time
t_mean = np.mean(t_lost)

Now, we can define the function $g(\beta)$ .

In [6]:

def root_fun(beta):
    exp_val = np.exp(beta * delta_t_ij)
    ret_val = -m * t_end - n * t_mean
    ret_val += np.sum(delta_t_ij * exp_val / (exp_val - 1))
    
    return ret_val

Now, we are left to find the MLE using SciPy's nifty solver scipy.optimize.fsolve(). We need to pass in the function whose root we are trying to find and also our initial guess at a root. We already have an approximate value for $\beta_\mathrm{MLE}$ , so we will use that.

In [7]:

beta_mle = scipy.optimize.fsolve(root_fun, beta_mle_approx)[0]

print(beta_mle)

0.013937036406431463

It's exactly the same as our approximate MLE, varying at the fourth decimal point. So, we are not too far off by using the approximate value for the MLE.

Checking our model¶

Now to check the model, we will plot the ECDF of the measured data along with the CDF of the theoretical Exponential distribution.

In [8]:

p = iqplot.ecdf(
    np.concatenate((t_lost, np.nan * np.ones(m))),
    q="time of bead loss (min)",
    y_range=[-0.025, 1.025],
)

# Theoretical curve
t_theor = np.linspace(0, 70, 200)
cdf = 1 - np.exp(-beta_mle * t_theor)
p.line(t_theor, cdf, line_width=2, line_color='orange')

bokeh.io.show(p)

We see that the agreement is quite poor. This could be due to a few reasons. The chemical process of cutting beads may not actually be well described by a Poisson process. This is very interesting! It would mean that there may be some more complex chemical dynamics at play, which we might want to explore with some more sophisticated experimentation.

Another possibility is that we made mistakes in determining how many beads were uncleaved. We filtered out stuck beads in our image processing protocol, but it is possible that there are beads with multiple DNA tethers, or those that may be stuck but moving on the surface.

If we lack confidence in our ability to detect tethered, but uncut beads, we could try another model, still assuming cutting is a Poisson process, but only taking into account cleaved beads.

An alternative model¶

The probability density function of bead loss times that we can observe is

\begin{aligned} (26) & f (t; β, t_{e n d}) \propto {\begin{array}{cll} β e^{- β t} & 0 \leq t \leq t_{e n d} \\ 0 & t > t_{e n d} . \end{array} \end{aligned}

$\begin{align} f(t ; \beta, t_\mathrm{end}) \propto \left\{\begin{array}{cll} \displaystyle{\beta\,\mathrm{e}^{-\beta t}} && 0 \le t \le t_\mathrm{end} \\[1em] 0 && t > t_\mathrm{end}. \end{array}\right. \end{align}$

We need to normalize this distribution. Let $\alpha$ be the constant or proportionality. Then,

\begin{aligned} (27) & \int_{0}^{\infty} d t f (t; β, t_{e n d}) = α \int_{0}^{t_{e n d}} d t β e^{- β t} = α (1 - e^{- β t_{e n d}}) = 1. \end{aligned}

$\begin{align} \int_0^\infty \mathrm{d}t\, f(t; \beta, t_\mathrm{end}) = \alpha\int_0^{t_\mathrm{end}} \mathrm{d}t\,\beta \mathrm{e}^{-\beta t} = \alpha\left(1 - \mathrm{e}^{-\beta t_\mathrm{end}}\right) = 1. \end{align}$

This gives $\alpha = \left(1 - \mathrm{e}^{-\beta t_\mathrm{end}}\right)^{-1}$ , and our probability distribution of cleavage times is

\begin{aligned} (28) & f (t; β, t_{e n d}) = {\begin{array}{cll} \frac{β e^{- β t}}{(1 - e^{- β t_{e n d}})} & 0 \leq t \leq t_{e n d} \\ 0 & t > t_{e n d} . \end{array} \end{aligned}

$\begin{align} &f(t; \beta, t_\mathrm{end}) = \left\{\begin{array}{cll} \displaystyle{\frac{\beta \mathrm{e}^{-\beta t}}{(1-\mathrm{e}^{-\beta t_\mathrm{end}})}} && 0 \le t \le t_\mathrm{end} \\[1em] 0 && t > t_\mathrm{end}. \end{array}\right. \end{align}$

So, for our set of $n$ observations, $\mathbf{t} \equiv \{t_1, t_2, \ldots t_n\}$ , we have a probability density function of

\begin{aligned} (29) & f (t; β, t_{e n d}) = \prod_{j = 1}^{n} f (t_{j}; β, t_{e n d}) . \end{aligned}

$\begin{align} f(\mathbf{t}; \beta, t_\mathrm{end}) = \prod_{j=1}^n f(t_j ; \beta, t_\mathrm{end}). \end{align}$

Writing this out, we have

\begin{aligned} (30) & f (t; β, t_{e n d}) = β^{n} \frac{e^{- n β \bar{t}}}{(1 - e^{- β t_{e n d}})^{n}}, \end{aligned}

$\begin{align} f(\mathbf{t} ; \beta, t_\mathrm{end}) = \beta^n\frac{\mathrm{e}^{-n\beta \bar{t}}}{(1-\mathrm{e}^{-\beta t_\mathrm{end}})^n}, \end{align}$

since all observed bead loss times are less than $t_\mathrm{end}$ . We will directly use this function to find the MLE for $\beta$ , as we showed in the previous model that the error in not explicitly taking into account the time interval of observation is small.

MLE for the alternative model¶

We wish to find the value of the parameter $\beta$ that leads to maximal $f(\mathbf{t};\beta, t_\mathrm{end})$ . We proceed as before by computing $\ln f(\mathbf{t};\beta, t_\mathrm{end})$ , differentiating, and setting to zero.

\begin{aligned} (31) & \ln f (t; β, t_{e n d}) = n \ln β - n β \bar{t} - n \ln (1 - e^{- β t_{e n d}}) . \end{aligned}

$\begin{align} \ln f(\mathbf{t};\beta, t_\mathrm{end}) = n\ln \beta - n\beta \bar{t} - n \ln\left(1 - \mathrm{e}^{-\beta t_\mathrm{end}}\right). \end{align}$

Differentiating with respect to $\beta$ gives

\begin{aligned} (32) & \frac{d}{d β} \ln f (t; β, t_{e n d}) = \frac{n}{β} - n \bar{t} - \frac{n t_{e n d} e^{- β t_{e n d}}}{1 - e^{- β t_{e n d}}} . \end{aligned}

$\begin{align} \frac{\mathrm{d}}{\mathrm{d}\beta}\,\ln f(\mathbf{t};\beta, t_\mathrm{end}) = \frac{n}{\beta} - n\bar{t} - \frac{nt_\mathrm{end}\mathrm{e}^{-\beta t_\mathrm{end}}}{1 - \mathrm{e}^{-\beta t_\mathrm{end}}} . \end{align}$

So, we are left to find the value of $\beta$ for which

\begin{aligned} (33) & g (β) = \frac{1}{β} - \bar{t} - \frac{t_{e n d} e^{- β t_{e n d}}}{1 - e^{- β t_{e n d}}} = 0. \end{aligned}

$\begin{align} g(\beta) = \frac{1}{\beta} - \bar{t} - \frac{t_\mathrm{end}\mathrm{e}^{-\beta t_\mathrm{end}}}{1 - \mathrm{e}^{-\beta t_\mathrm{end}}} = 0. \end{align}$

We can again use scipy.optimize.fsolve() to find the value of the MLE. We will use $1/\bar{t}$ as our initial guess.

In [9]:

def alt_root_fun(beta):
    exp_val = np.exp(-beta * t_end)

    return 1 / beta - t_mean - t_end * exp_val / (1 - exp_val)


# Find the MLE!
alt_beta_mle = scipy.optimize.fsolve(alt_root_fun, 1 / t_mean)[0]

# Report result
print(alt_beta_mle)

0.041951259352945

We can also get the time scale.

In [10]:

print(1 / alt_beta_mle)

23.837186664333107

This is much shorter, about 24 minutes. This might be expected, as we could have over-assigned beads that we did not lose.

Plotting the CDF¶

Now that we have a maximum likelihood estimate for $\beta$ , we can compute the theoretical CDF and see how it looks compared to our measured ECDF. The theoretical CDF is

\begin{aligned} (34) & CDF = \int_{0}^{t} d t^{'} f (t^{'}; β, t_{e n d}) = \frac{1 - e^{- β t}}{1 - e^{- β t_{e n d}}} . \end{aligned}

$\begin{align} \text{CDF} = \int_0^t\mathrm{d}t'\,f(t'; \beta, t_\mathrm{end}) = \frac{1-\mathrm{e}^{-\beta t}}{1 - \mathrm{e}^{-\beta t_\mathrm{end}}}. \end{align}$

We can code this up and plot it along with our empirical CDF.

In [11]:

p = iqplot.ecdf(t_lost, q="time of bead loss (min)")

# Compute theoretical CDF
t_theor = np.linspace(0, t_end, 400)
cdf = (1 - np.exp(-alt_beta_mle * t_theor)) / (1 - np.exp(-alt_beta_mle * t_end))
p.line(t_theor, cdf, line_width=2, line_color='orange')

bokeh.io.show(p)

This looks much better, suggesting that perhaps we did misidentify beads we did not lose. It's a pity, since that limits the data we can use. This motivates us to collect images as long as possible.