The Lincoln Index Problem¶

License: Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0)

# If we're running on Colab, install libraries

import sys
IN_COLAB = 'google.colab' in sys.modules

if IN_COLAB:
    !pip install pymc3
    !pip install arviz
    !pip install graphviz

In an excellent blog post, John D. Cook wrote about the Lincoln index, which is a way to estimate the number of errors in a document (or program) by comparing results from two independent testers.

Here's his presentation of the problem:

"Suppose you have a tester who finds 20 bugs in your program. You want to estimate how many bugs are really in the program. You know there are at least 20 bugs, and if you have supreme confidence in your tester, you may suppose there are around 20 bugs. But maybe your tester isn't very good. Maybe there are hundreds of bugs. How can you have any idea how many bugs there are? There's no way to know with one tester. But if you have two testers, you can get a good idea, even if you don't know how skilled the testers are."

Suppose the first tester finds 20 bugs, the second finds 15, and they find 3 in common; how can we estimate the number of bugs?

I'll use the following notation for the data:

For this model, I'll express the data in different notation:

k11 is the number of bugs found by both testers,
k10 is the number of bugs found by the first tester but not the second,
k01 is the number of bugs found by the second tester but not the first, and
k00 is the unknown number of undiscovered bugs.

Here are the values for all but k00:

k10 = 20 - 3
k01 = 15 - 3
k11 = 3

num_seen = k01 + k10 + k11
num_seen

32

Here's the model:

import pymc3 as pm

with pm.Model() as model5:
    p0 = pm.Beta('p0', alpha=1, beta=1)
    p1 = pm.Beta('p1', alpha=1, beta=1)
    N = pm.DiscreteUniform('N', num_seen, 350)
    
    q0 = 1-p0
    q1 = 1-p1
    ps = [q0*q1, q0*p1, p0*q1, p0*p1]
    
    k00 = N - num_seen
    data = pm.math.stack((k00, k01, k10, k11))
    y = pm.Multinomial('y', n=N, p=ps, observed=data)

with model5:
    trace5 = pm.sample(500)

Multiprocess sampling (2 chains in 2 jobs)
CompoundStep
>NUTS: [p1, p0]
>Metropolis: [N]
Sampling 2 chains, 0 divergences: 100%|██████████| 2000/2000 [00:01<00:00, 1577.49draws/s]
The rhat statistic is larger than 1.05 for some parameters. This indicates slight problems during sampling.
The estimated number of effective samples is smaller than 200 for some parameters.

with model5:
    pm.plot_posterior(trace5)

with model5:
    pm.traceplot(trace5)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-6-5ba0275a3a86> in <module>
      1 with model5:
----> 2     pm.traceplot(trace5)

~/anaconda3/envs/ThinkBayes2/lib/python3.8/site-packages/pymc3/plots/__init__.py in wrapped(*args, **kwargs)
     21                 warnings.warn('Keyword argument `{old}` renamed to `{new}`, and will be removed in pymc3 3.8'.format(old=old, new=new))
     22                 kwargs[new] = kwargs.pop(old)
---> 23             return func(*args, **kwargs)
     24     return wrapped
     25 

~/anaconda3/envs/ThinkBayes2/lib/python3.8/site-packages/pymc3/plots/__init__.py in traceplot(*args, **kwargs)
     39     try:
     40         kwargs.setdefault('compact', True)
---> 41         return az.plot_trace(*args, **kwargs)
     42     except TypeError:
     43         kwargs.pop('compact')

~/anaconda3/envs/ThinkBayes2/lib/python3.8/site-packages/arviz/plots/traceplot.py in plot_trace(data, var_names, coords, divergences, figsize, rug, lines, compact, combined, legend, plot_kwargs, fill_kwargs, rug_kwargs, hist_kwargs, trace_kwargs, backend, backend_kwargs, show)
    216 
    217     plot = get_plotting_function("plot_trace", "traceplot", backend)
--> 218     axes = plot(**trace_plot_args)
    219 
    220     return axes

~/anaconda3/envs/ThinkBayes2/lib/python3.8/site-packages/arviz/plots/backends/matplotlib/traceplot.py in plot_trace(data, var_names, divergences, figsize, rug, lines, combined, legend, plot_kwargs, fill_kwargs, rug_kwargs, hist_kwargs, trace_kwargs, plotters, divergence_data, colors, backend_kwargs, show)
    129 
    130         if len(value.shape) == 2:
--> 131             _plot_chains_mpl(
    132                 axes,
    133                 idx,

~/anaconda3/envs/ThinkBayes2/lib/python3.8/site-packages/arviz/plots/backends/matplotlib/traceplot.py in _plot_chains_mpl(axes, idx, value, data, colors, combined, xt_labelsize, rug, trace_kwargs, hist_kwargs, plot_kwargs, fill_kwargs, rug_kwargs)
    262         if not combined:
    263             plot_kwargs["color"] = colors[chain_idx]
--> 264             plot_dist(
    265                 values=row,
    266                 textsize=xt_labelsize,

~/anaconda3/envs/ThinkBayes2/lib/python3.8/site-packages/arviz/plots/distplot.py in plot_dist(values, values2, color, kind, cumulative, label, rotated, rug, bw, quantiles, contour, fill_last, textsize, plot_kwargs, fill_kwargs, rug_kwargs, contour_kwargs, contourf_kwargs, pcolormesh_kwargs, hist_kwargs, ax, backend, backend_kwargs, show, **kwargs)
    151             hist_kwargs = {}
    152 
--> 153         hist_kwargs.setdefault("bins", get_bins(values))
    154         hist_kwargs.setdefault("cumulative", cumulative)
    155         hist_kwargs.setdefault("color", color)

~/anaconda3/envs/ThinkBayes2/lib/python3.8/site-packages/arviz/plots/plot_utils.py in get_bins(values)
    129     bins_fd = 2 * iqr * values.size ** (-1 / 3)
    130 
--> 131     width = round(np.max([1, bins_sturges, bins_fd])).astype(int)
    132 
    133     return np.arange(x_min, x_max + width + 1, width)

AttributeError: 'int' object has no attribute 'astype'