核酸检测多少人为一组混检合适?

今天在想一个有趣的问题。核酸检测的混检,必然是当患病率越低时,则混到一个管子的人数越多越好,因为这样检测的期望次数就会更少。那么问题来了,当患病率多高时,混检就失去了意义?当混检没失去意义时,检测多少次才会令期望的检测次数最低。

知乎上有人给出了这个问题的 答案。然而,这个答案我并不是非常地满意,他们只是考虑总人数比较多的时候,用样本概率去近似估计条件概率,有一定的近似成分在里面。能精确逼近的,为什么要用估计?我不喜欢这个。某种意义上来说,他们给出的答案是 “不对” 的。“不对” 加引号是因为,即使不对,在样本数足够大的情况下,几乎是对的。

来,我们来考虑一下这个问题。不妨假设总人数为 N N N,患病(敏感词汇,用患病替代)的人数为 M M M,假设混检的策略采用 k k k 个人一组进行混检,检出患病,这 k k k 个人则进行一次额外的检测。问,给定 N N N M M M 的情况下,总的检测次数(期望意义下)和 k k k 是什么样的一个关系?为了方便,假设下述的所有相除数量关系都是可以整除的,这个不是重要的。

如果我们以所有人数作为考虑,计算他们总检测次数的期望,你会发现,它的检测次数出现的情况非常多,小到 N / k + k N/k+k N/k+k,大到 N / k + M k N/k+Mk N/k+Mk。而要计算每个次数出现的概率,要算很多组合数,非常麻烦。因为,我们不妨换一个角度,考虑每个人消耗检测次数的期望值,再把它乘以总人数 N N N,那么就得到这个整体的检测次数的期望。

对于一个人来说,他消耗的检测次数要么是 1 / k 1/k 1/k,要么是 1 + 1 / k 1+1/k 1+1/k,只有这两种情况,相对简单多了。 1 / k 1/k 1/k 当且仅当他的这一组里面,没有人患病。它的概率是:
P ( 1 / k ) = C N − M k C N k = ∏ i = 0 k − 1 N − M − i N − i P(1/k) = \frac{C_{N-M}^k}{C_N^k}=\prod_{i=0}^{k-1}\frac{N-M-i}{N-i} P(1/k)=CNkCNMk=i=0k1NiNMi
同理,只要这一管子里面有人患病,那么他消耗的检测次数就会变成 1 + 1 / k 1+1/k 1+1/k。这种情况发生的概率为:
P ( 1 + 1 / k ) = 1 − C N − M k C N k = 1 − ∏ i = 0 k − 1 N − M − i N − i . P(1+1/k) = 1- \frac{C_{N-M}^k}{C_N^k}=1-\prod_{i=0}^{k-1}\frac{N-M-i}{N-i}. P(1+1/k)=1CNkCNMk=1i=0k1NiNMi.
那么,一个人的消耗的检测次数的数学期望为:
E ( k ) = 1 / k ∗ P ( 1 / k ) + ( 1 + 1 / k ) ∗ P ( 1 + 1 / k ) = k + 1 k + ∏ i = 0 k − 1 M − N + i N − i E(k) = 1/k*P(1/k)+(1+1/k) *P(1+1/k) =\frac{k+1}{k}+\prod_{i=0}^{k-1}\frac{M-N+i}{N-i} E(k)=1/kP(1/k)+(1+1/k)P(1+1/k)=kk+1+i=0k1NiMN+i
MMP,这个似乎没办法再化简。不是我不行,我用 MATLAB 符号计算试了一下,确实没法化简。用到的代码如下:

clc
clear
close all
syms N M 
k = 5;
Pk1 = 1;
for i = 0:k-1Pk1 = Pk1*(N-M-i)/(N-i);
end
Pk1_ = 1-Pk1;
E = 1/k*Pk1+(1+1/k)*Pk1_;
Esim = simplify(E)
Efac = factor(Esim);
latex(Esim)
latex(Efac)

那么,总的检测次数的数学期望就是 N ∗ E ( k ) N*E(k) NE(k),一般我们考虑 E ( k ) E(k) E(k) 就可以了。我们用两种方法来验证一下这个结果,第一个方法是找一些简单的值和手算做对比,第二个方法是用程序做一个大数据量下的仿真。

M = 1 M=1 M=1 N = 4 N=4 N=4 k = 2 k=2 k=2 时, N ∗ E ( k ) = 4 N*E(k)=4 NE(k)=4,和我们手算的结果一样。所用到的程序如下。

NE(M,N) = N*E;
NEfun = matlabFunction(NE);
NEfun(2,6)

M = 2 M=2 M=2 N = 6 N=6 N=6 k = 3 k=3 k=3 时, N ∗ E ( k ) = 6.8 N*E(k)=6.8 NE(k)=6.8,和我们手算的结果 6.5 不一样。手算的过程为, 6 6 6 个人两组,至少要做两次。两个患病的进入一个管子的概率如同时抛两枚硬币,硬币朝向是相同的概率为 0.5。故而,有 0.5 的可能要多做 3 次,有另外 0.5 的可能要多做 6 次。一次总的次数 2 + 3 ∗ 0.5 + 6 ∗ 0.5 = 6.5 2+3*0.5+6*0.5 = 6.5 2+30.5+60.5=6.5

这里的结果发生了不一致,到底是谁对谁错呢?让我们用程序仿真来模拟一下。

clc
clear
close all
M = 2;
N = 6;
k = 3;
Ms = 1:M;s = 0;
NN = 1000000;
N_k = N/k;
for ii=1:NNR = reshape(randperm(N),k,[]);res = 0;for i=1:N_kif(intersect(R(:,i),Ms)>0)res = res+1;endends = s+res*k+N_k;s/ii
end

模拟的结果是 6.8,说明我们手算的那个思路是不对的,我来看看为什么。其实很简单,2个管子,6 个人,每个管子 3 个人,2 个患病的人进入同一个管子的概率不是 1/2,而是 1 / 2 ∗ 2 / 5 ∗ 2 = 2 / 5 1/2*2/5*2=2/5 1/22/52=2/5。这个和抛硬币是不一样的,仔细想想就知道了,这里不再赘述。一言以蔽之,抛硬币没有一管三人的约束,是不一样的。

接下来,我们来看一下,在 M M M N N N 固定的情况下, E E E k k k 是什么样一个限制关系,以及 E ( k ) E(k) E(k) 的极值点在哪。上面的论述都是基于 k > 1 k>1 k>1 的时候,当 k = 1 k=1 k=1 的时候,因为就一个人,即使测出患病,也没必要再测一次。所以,需要对公式做一个修正。考虑
E ( k ) = k + 1 k + ∏ i = 0 k − 1 M − N + i N − i , k > 1 E ( k ) = 1 , k = 1. E(k) =\frac{k+1}{k}+\prod_{i=0}^{k-1}\frac{M-N+i}{N-i},k>1\\ E(k) =1,k=1. E(k)=kk+1+i=0k1NiMN+ik>1E(k)=1k=1.

事实上, E E E 关于 k k k 求导 E ′ ( k ) E'(k) E(k) 是一件非常麻烦的事情,不是初等数学可以做的。所以,我们寄希望于用 MATLAB 画一下图,再看看结果。为了把细节看清楚,我分了两张图画。
在这里插入图片描述
从图上可以看得出来,当患病率高于 0.3 的时候,就已经不适合混检了。
请添加图片描述
可以看得出来,患病率越低,可以达到的极值点越小。且随着 k k k 的增加, E ( k ) E(k) E(k) 先减后增,最后总是趋向于 1 的。当患病率低到 0.01 左右的时候,差不多 10 个人混检是比较合理的。

clc;clear;close all;
Rs = [0.05 0.04 0.02 0.011 0.01 0.005 0.001 0.0005];
%Rs = 4/100000;
for R = RsN = 40000000000000;M = round(N*R);K = 50;for k=1:Ky(k) = E(N,M,k);endplot(1:K,y);title('个人消耗检测数随混检人数 k 的变化')xlabel('混检人数 k')ylabel('每个人消耗试剂数期望 E(k)')hold on;[val,ind] = min(y);text_points(ind,val,val)%scatter(ind,val)
end
for i = 1:length(Rs)legends{i} = ['患病率R=' num2str(Rs(i))];
end
legend(legends)function Ev = E(N,M,k)
if(k==1)Ev=1;return;
end
Pk1 = 1;
for i = 0:k-1Pk1 = Pk1*(N-M-i)./(N-i);
end
Pk1_ = 1-Pk1;
Ev = 1./k*Pk1+(1+1./k).*Pk1_;
endfunction text_points(x,y,R)
sz = length(x);
if(size(x,2)>1)x = x.';
end
if(size(y,2)>1)y = y.';
end
if(size(R,2)>1)R = R.';
end
strs = [num2str(x) ',' num2str(R)];
text(x,y,strs);
end

还有一个问题是,考虑到北京市当前的极低的患病率,大约 4/100000,那么得到的最优的 k k k 就应该是 160。但是现在却不以这个数值执行,这是为什么呢?这是因为取这个最优值消耗的试剂虽然是少了。但是对于被检测者来说,不管是混检还是单检,做一次核酸真的就是要付出做一次核酸的成本,排一次队,刮一次鼻子,遭一次罪。混检对测试成本消耗来说,确实是 1 / k 1/k 1/k,但对于被测试这来说,他就是 1。极端情况,如果 k = 1 k=1 k=1 ,百姓只要做一次核酸,核酸免费的时代,他们当然更愿意单检。反之,如果 k k k 太大了,虽然检测次数期望达到了最优,但是一旦某个管子检测出患病,哪怕里面只有一个人患病,就要白白拉着其他 k − 1 k-1 k1 个人再做一次检测。将被检测者的感受纳入考量,本着以人为本,所以这个 k k k 不一定要取到最优,有时候牺牲一点测试成本,换来大家的轻松,也是很好的选择。

有人说,可以再采用的时候,把每个人的样本分两份,当检测出患病的时候,再把备份拿出来单独测,这样就不需要喊人来车第二次了。考虑到患病率极低,这种分两份的人工成本是极高的,已经远远超过了把一起混检的人喊来再测一次的成本,更远远超过了不采用最优参数 k k k 而采用更小的,所损失的测试成本。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/20557.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

用python写了个全国疫情中高风险地区查询

最近用python写了个全国疫情中高风险地区查询的爬虫代码,分享给大家一起交流,希望得到不同思路的指教,让代码更简洁,运行效率更高。 总体思路 1、找到可供查询的源网站 2、分析、获取查询的API 3、构造API 4、获取全国中高风险地…

如何查询澳大利亚药监局(TGA)药品信息数据

澳大利亚医疗用品管理局(Therapeutic Goods Administration),简称TGA,是澳洲医疗用品的监管机构,负责一系列评估和监管确保澳洲药品保质保量。TGA监管的产品范围包括药品,医疗器械,血液及血液产品等。 澳大利亚药品分…

文字识别核酸检测结果并导出Excel

python下文字识别核酸检测报告信息并生成Excel数据表 前言:疫情下,隔几天就需要做一次核酸检测,核酸检测截图的收取工作称为部门工作的难点,参考最近热搜话题复旦大学辅导员用python搞定核酸检测核查难的问题,也自己动…

新库上线 | CnOpenData中国核酸检测机构及采样点数据

中国核酸检测机构及采样点数据 一、数据简介 2020年1月21日,国家卫健委发布1号公告,将新型冠状病毒感染的肺炎纳入《中华人民共和国传染病防治法》规定的乙类传染病,并采取甲类传染病的预防、控制措施。目前,新型冠状病毒肺炎防控…

电脑蓝牙与蓝牙适配器使用

这篇文章一开始写的很水,但最近重装电脑还是遇到了这个问题。下面更新一下,解决搜不到的问题。 ------------------------------------------------------------------------------------------------------ 1.确保你的手机能够蓝牙连接到耳机&#xf…

chatgpt赋能python:Python设置画布背景图——让你的图像更具美感

Python设置画布背景图——让你的图像更具美感 Python是一门流行的编程语言,被广泛使用于数据分析、科学计算和图像处理等多个领域。在图像处理方面,Python使用matplotlib作为主要的作图工具,matplotlib带来了许多方便易用的工具,…

为什么ChatGPT用强化学习而非监督学习?

为什么ChatGPT非得用强化学习,而不直接用监督学习?原因不是那么显而易见。在上周发布的《John Schulman:通往TruthGPT之路》一文中,OpenAI联合创始人、ChatGPT主要负责人John Schulman分享了OpenAI在人类反馈的强化学习&#xff0…

适配PyTorch FX,OneFlow让量化感知训练更简单

作者 | 刘耀辉 审稿 | BBuf、许啸宇 1 背景 近年来,量化感知训练是一个较为热点的问题,可以大大优化量化后训练造成精度损失的问题,使得训练过程更加高效。 Torch.fx在这一问题上走在了前列,使用纯Python语言实现了对于Torch.nn.M…

mongodb charts对mongodb数据进行分析和展示

mongodb charts 安装教程 安装环境什么是mongodb charts下载mongodb charts等准备工作配置mongodb charts创建用户启动和停止MongoDB图表故障排除web展示 安装环境 系统环境:ubuntu 16.04 docker 版本:Docker version 18.09.0 mongo 版本:Mo…

小白量化彩票实战(4)彩票特征号码重号、邻号、连号和表格展示

小白量化彩票实战(4)彩票特征号码重号、邻号、连号和表格展示 我写彩票的博客,不是鼓励大家去买彩票,读者要以学习编程和娱乐的思想来看待。兴趣是学习最大的动力! 彩票的号码特征很多,我们本篇介绍几个简单的号码特征…

一个小把戏算法,获取大乐透,并且计算出最佳的结果(Qt C++ 和Android共用)

无聊的国庆,总得做点什么好玩的是不是,那就写代码获取大乐透,让后按照自己的算法推测下一期的结果吧。 话不多说,上代码 Widget::Widget(QWidget *parent): QWidget(parent), ui(new Ui::Widget) {ui->setupUi(this);initDat…

算力军备竞赛白热化 “卖铲人”联想集团竞争力如何?

继微软通过OpenAI推出GPT系列、谷歌推出Bard和PaLM-E2之后,国内AI大模型也呈百家争鸣态势,年初至今,国内科技巨头几乎都发布了自研AI大模型产品,AI竞赛全面升级的背后,是全球科技巨头们对算力的争夺,作为算…

chatgpt赋能python:Python制图中如何添加文字

Python 制图中如何添加文字 介绍 制图通常不仅需要展示数据,还需要向读者传递信息。而添加文字是一种直接有效的方式,可以帮助读者更好地理解图表。 Python 图形库众多,如 Matplotlib、Seaborn、Plotly 等,它们都提供了向图表中…

【送书福利-第八期】《硅基物语.AI大爆炸: ChatGPT→AIGC→GPT-X→AGI进化→魔法时代→人类未来》

大家好,我是洲洲,欢迎关注,一个爱听周杰伦的程序员。关注公众号【程序员洲洲】即可获得10G学习资料、面试笔记、大厂独家学习体系路线等…还可以加入技术交流群欢迎大家在CSDN后台私信我! 本文目录 一、前言二、内容介绍三、作者介…

算法工程师体验了一下chatGPT,已经上瘾了!

chatGPT持续刷屏,作为能写代码,能修bug的超级工具,CV君必须体验一把! 首先来一个基本操作,让chatGPT写一段Python程序,使用YOLOv5对图像中的目标进行检测,找出有狗没有猫的图片: 对YOLOv5这种公…

还有人不懂 ChatGPT,不焦虑吗?(文末赠书)

,不 如果有一本书 可以让人理解“AI大爆炸”新纪元 那就是《碳基物语》 半年以来,ChatGPT点燃文明新火把 对AIGC和AGI的讨论也甚嚣尘上‍‍‍‍ AI会取代人类吗? 人工智能会拥有智慧吗? ChatGPT到底该怎么玩? 我该如何…

一想到还有95%的人不懂ChatGPT,我就焦虑了

如果有一本书 可以让人理解“AI大爆炸”新纪元 那就是《碳基物语》 半年以来,ChatGPT点燃文明新火把 对AIGC和AGI的讨论也甚嚣尘上‍‍‍‍ AI会取代人类吗? 人工智能会拥有智慧吗? ChatGPT到底该怎么玩? 我该如何利用AIGC提升生产…

5月书讯 | 《这就是ChatGPT》来了!

叮~又到了书讯时间,本月好书众多,姗姗来迟。 在这个数字化的时代,我们仍然相信纸质书的魅力,可以让人沉静下来,回归到阅读的本质。五月盛夏伊始,炎炎夏日,我们精心挑选了 10 本好书,…

通过AI的自白,开启ChatGPT学习之旅!

如果有一本书 可以让人理解“AI大爆炸”新纪元 那就是《碳基物语》 半年以来,ChatGPT点燃文明新火把 对AIGC和AGI的讨论也甚嚣尘上‍‍‍‍ AI会取代人类吗? 人工智能会拥有智慧吗? ChatGPT到底该怎么玩? 我该如何利用AIGC提升生产…

Anaconda安装的python环境中“No module named pip” 和 “ ‘pip‘ is a package and cannot be directly executed”问题

一. 没有pip3问题 找到安装anaconda的文件夹,点击Scripts(利用anaconda安装的python虚拟环境都在这里),确定是否存在一个easy_install.exe的程序,如果有请往下看,如果没有进入直接进入第4步。 打开 Anaconda Prompt 或 cmd &…