新冠病毒的组合检测

近日，武汉市组织全市一千万余人进行新冠病毒核酸检测。很多认为对一千万人检测就必然用掉一千万份试剂，其实不然：在阳性率不高的情况下，运用组合检测的方法能大大降低检测所需要的试剂量。具体来说，我们在检测的时候可以把多人的样本融合在一起，然后检测整组是否为阳性。如整组为阴性，那么所有人应当都为阴性；如整组为阳性，那么再对组内的每一个人进行检测。

我们假设人口为$n$，阳性率为$a$，在分组时每$k$人一组。现有如下事实：

  • 共有$\left\lceil\frac{n}{k}\right\rceil=m$组
  • 共有$na$人阳性，$n(1-a)$人阴性

假设所有人被随机分配到了这$m$组内，对于组内的每一个人，其为阴性的概率是$(1-a)$。显然地，某一组至少存在一个阳性的概率为$1-(1-a)^k$。对于每一个组，检测量的数学期望为：

\begin{align*} (1-a)^k \cdot 1 + \left[ 1- (1-a)^k \right] \cdot k \end{align*}

由于一共有$m$个组，最后总共的检测量为：

\begin{align*} m\left[ (1-a)^k + k - k(1-a)^k \right] \end{align*}
In [5]:
import numpy as np
import matplotlib
matplotlib.rcParams['figure.dpi'] = 200
matplotlib.rcParams['font.sans-serif'] = ['FandolSong']  # 用来正常显示中文标签
matplotlib.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
import matplotlib.pyplot as plt


def get_num_total_test(n, a, k):
    num_group = np.ceil(n / k) #m
    neg_power = np.power(1.0-a, k)
    test_per_group = neg_power + (1 - neg_power) * k
    return num_group * test_per_group
In [2]:
get_num_total_test(10000000, 0.05, 64)
Out[2]:
9630621.754640577
In [9]:
n_ = 10000000
a_s_ = [0.5, 1, 2, 3, 5, 10] # in percentage
k_s_ = [2, 4, 8, 16, 32, 64]

def graph_total_tests():
    global n_, a_s_, k_s_
    n = n_
    a_s = np.asarray(a_s_, np.float) / 100.0
    k_s = np.asarray(k_s_)
    
    bar_total_width = 0.8
    num_bar_per_index = a_s.size
    bar_width = bar_total_width / num_bar_per_index
    
    # compute number of tests, grouped by each a
    num_tests = []
    for a in a_s:
        tests = get_num_total_test(n, a, k_s)
        num_tests.append(tests)
    
    fig, ax = plt.subplots(1, 1)
    for i in range(len(num_tests)):
        cur_data = num_tests[i]
        ax.bar(np.arange(k_s.size) - bar_total_width / 2 + i * bar_width, cur_data, bar_width, 
               label='$a={:.1f}\%$'.format(a_s[i] * 100.0))
    
    ax.set_xlabel('每组人数（$k$）')
    ax.set_xticks(np.arange(k_s.size))
    ax.set_xticklabels(k_s)
    ax.set_ylabel('总检测量')
    ax.ticklabel_format(style='plain', axis='y')
    fig.suptitle('分组检测参数与总检测量（$n=10^7$）')
    
    ax.legend(ncol=2, title='阳性率')
    
graph_total_tests()
In [ ]: