基于Python实现的分组检测
引言
分组检测是当今医学中经常出现的一种对疾病检测的有效方法。对于收集咽拭子样本并进行检测这个操作是需要一定成本的,尤其是疫情突然爆发,在短时间内,我们需要在当地众多居民中检测出携带病毒的个体,在新冠患病率不太高的地区,如果我们对众多个体进行逐个检测,可定会造成很大的浪费,因为其中大部分人是不需要检测的。对于这个问题,分组检测可以有效解决。运用python编程技术和概率论相关知识可以轻松验证数学期望的合理性
1.问题提出
在一个人数很多的团体中检测新型冠状病毒的携带者,有N个人进行测试,该集体中患病率为p,每个样本是否患病是相互独立的,现有两种检测方法:(1)逐个检测(2)按x个人一组进行分组,把x个人的咽拭子样本混合在一起,如果混合样本呈现阳性,就说明在这x个人中有人患病,然后把这x个人的样本再检测一次。如果x个人中没有人患病,则该组只需检测1次,若有人患病,则需检测x+1次。
(1)验证化验次数的数学期望函数的合理性
(2)分组人数取何值(在某有效范围内)时,节省的工作效率最高
(3)当患病率控制在什么范围内,分组检测的方法有效
2.验证化验次数数学期望函数的合理性
2.1样本来源:
2020年10月11日,青岛市发现3例新冠肺炎无症状感染者,随后组织开展大规模流调排查和分类检测。截至2020年11月16日18时,共收集样本10899145份。
假设这108万样本分5轮进行检测,有10个检测点,每次检测21800份,青岛地区患病率为0.01。
2.2方法简述
本次取证通过对比两条曲线(x为分组人数,y为检测次数)作比较得出结论。曲线一:以2-100的连续整数作为x,通过期望计算公式计算y。曲线二:以2-100之间的连续整数作为x,并根据x进行分组:如果21800%x=0,则共有21800/x组。如果21800%x!=0,则共有21800/x+1组,最后一组的人数为21800%x。利用random包生成与权重有关的随机样本(在样本中,测试结果为阳性的是1,结果为阴性的是1)。如果在一组中包含1,则该组的检测次数为x+1,反之,检测次数为1。每组检测次数相加即为x对应的y.
2.3代码展示
```python
def count_exception(sumSamples,mobidity,x): #计算期望的函数
f=(sumSamples//x)*(1*(1-mobidity)**x+(x+1)*(1-(1-mobidity)**x))
return f
'''青岛核酸检测'''
import random
import matplotlib.pyplot as plt
import numpy as np
sumSamples=21800
mobidity=0.01
def weighted_random(weights): #权重随机数
number=random.random()*sum(weights.values())
for x,p in weights.items():
if number<p:
break
return x
def creat_samples(num,weights): #根据权重产生样本
samples=[]
for i in range(num):
samples.append(weighted_random(weights))
return samples
weights={1:mobidity,0:1-mobidity}
samples=creat_samples(sumSamples, weights)
#print(samples)
#ro_visualization
'''
factors=[factor for factor in range(2,100) if not sumSamples%factor]
x=np.array(factors)
f=count_exception(sumSamples, mobidity, x)
print(x)
print(f)
print(type(x),type(f))
plt.plot(x, f,'ro')
'''
#产生分组人数为range(2,100)区间的公式化期望函数
factors1=[factor for factor in range(2,100)]
x=np.array(factors1)
f1=count_exception(sumSamples, mobidity, x)
plt.plot(f1, 'b')
#产生分组人数为range(2,100)区间的随机产生样本曲线
samples=creat_samples(sumSamples, weights)
#print(len(samples))
#print(samples)
#numGroup=sumSamples//factors1 if not sumSamples%factors1 else 'error'
def sum_test(samples,fi): #fi为factors中的一个元素numEachGroup
numGroup=sumSamples//fi
sumTest=0
minTest=1
maxTest=fi+1
#print(numGroup)
for group in range(numGroup):
startPos=group*fi
endPos=startPos+fi
eachGroupSample=samples[startPos:endPos]
#print('开始位置与结束位置:',startPos,endPos)
#print("每组中的元素值:",eachGroupSample)
if 1 in eachGroupSample:
sumTest+=maxTest
else:
sumTest+=minTest
if sumSamples%fi:
m=sumSamples%fi
eachGroupSample=samples[len(samples)-m:len(samples)]
if 1 in eachGroupSample:
sumTest+=sumSamples%fi+1
else:
sumTest+=minTest
return sumTest
sumTest=[]
factors=[factor for factor in range(2,100)]
for i in factors1:
sumTest.append(sum_test(samples, i))
#print(sumTest)
plt.plot(factors1, sumTest,'r')
plt.show()
2.4结果及分析
蓝色曲线为根据分组检测数学期望函数计算所得的分组人数与检测次数之间的关系;
红色权限为根据分组人数产生随机样本并计算检测次数得到的曲线关系
由图像可得:两条曲线非常相近,检验次数数学期望计算公式与现实规律基本相符,具有合理性。
3.探究分组人数上限及最优分组人数
3.1方法简介
已知样本总量,患病率、期望最大工作量(%)。分组人数为2-100之间的连续整数,利用期望公式计算检测次数的期望。期望与样本总数做比值即为实际工作量,如果该实际工作量大于期望最大工作量则返回最大分组人数。
3.2代码展示
def count_exception(sumSamples,mobidity,x):
f=(sumSamples//x)*(1*(1-mobidity)**x+(x+1)*(1-(1-mobidity)**x))
return f
#分组上限函数(max_factors) 注:rate参数
def max_numGroup(sumSamples,mobidity,rate):
factors=[factor for factor in range(2,100)] #在一组连续的数中
max_numEachGroup=0
for i in factors:
test_rate=count_exception(sumSamples, mobidity, i)//sumSamples
#print(test_rate)
if test_rate>rate:
max_numEachGroup=i
break
else:
max_numEachGroup=i
return max_numEachGroup
print(max_numGroup(21800,0.1,0.5))
mobidity=0.1
sumSamples=21800
import numpy as np
import matplotlib.pyplot as plt
factors=[factor for factor in range(2,max_numGroup(sumSamples, mobidity, 0.5)) if not sumSamples%factor]
x=np.array(factors)
f=count_exception(sumSamples, mobidity, x)
print(x)
print(f)
print(type(x),type(f))
plt.plot(x, f,'ro')
3.3结果及分析
当样本总人数为21800,患病率为0.1,期望最大工作量为0.5,得到分组人数上限为34.
分组人数为2-34之间且能被总样本人数整除的整数,分组人数的列表为[ 2 4 5 8 10 20 25 40 50],分组人数与检测次数之间的散点图。
当分组人数为4时,检测次数最少。
4.患病率控制在什么范围内,分组检测的方法有效
4.1方法简介
以患病率作为变量,分组人数为2-100之间的整数,运用检测次数数学期望函数计算检测次数。形成患病率为x变量,人均检测次数为y变量的曲线进行比较。人均检测次数为总检测次数与样本总量的比值(当人均检测次数大于1时,分组检测方法无效;当人均检测次数小于1时,分组检测次数有效)
4.2代码展示
import random
def count_exception(sumSamples,mobidity,x):
f=(sumSamples//x)*(1*(1-mobidity)**x+(x+1)*(1-(1-mobidity)**x))
return f
def weighted_random(weights):
number=random.random()*sum(weights.values())
for x,p in weights.items():
if number<p:
break
return x
def creat_samples(num,weights):
samples=[]
for i in range(num):
samples.append(weighted_random(weights))
return samples
def max_numGroup(sumSamples,mobidity,rate):
factors=[factor for factor in range(2,100)] #在一组连续的数中
max_numEachGroup=0
for i in factors:
test_rate=count_exception(sumSamples, mobidity, i)//sumSamples
#print(test_rate)
if test_rate>rate:
max_numEachGroup=i
break
else:
max_numEachGroup=i
return max_numEachGroup
import matplotlib.pyplot as plt
mobidity_list=[0.05,0.2,0.35,0.5,0.65,0.8,0.95]
sumSamples=10000
rate=0.7
def sum_test(samples,fi): #fi为factors中的一个元素numEachGroup
numGroup=sumSamples//fi
sumTest=0
minTest=1
maxTest=fi+1
#print(numGroup)
for group in range(numGroup):
startPos=group*fi
endPos=startPos+fi
eachGroupSample=samples[startPos:endPos]
#print('开始位置与结束位置:',startPos,endPos)
#print("每组中的元素值:",eachGroupSample)
if 1 in eachGroupSample:
sumTest+=maxTest
else:
sumTest+=minTest
if sumSamples%fi:
m=sumSamples%fi
eachGroupSample=samples[len(samples)-m:len(samples)]
if 1 in eachGroupSample:
sumTest+=sumSamples%fi+1
else:
sumTest+=minTest
return sumTest
sumTest=[]
factors1=[factor for factor in range(2,100) if not sumSamples%factor]
for i in factors1:
sumTest.append(sum_test(sumSamples, i))
#print(sumTest)
plt.plot(factors1, sumTest,'r')
plt.show()
for mobidity in mobidity_list:
saveRate_list=[]
save_rate=0
weights={1:mobidity,0:1-mobidity}
samples=creat_samples(sumSamples, weights)
factors=[i for i in range(2,max_numGroup(sumSamples, mobidity, rate)) if not sumSamples%i]
print(factors)
for fi in factors:
sumTest=sum_test(samples,fi)
save_rate=sumTest/sumSamples
saveRate_list.append(sumTest)
print(saveRate_list)
plt.plot(factors,saveRate_list)
4.3结果及分析
样本总量为10000,分组人数为2-100之间的整数
蓝橙绿红蓝曲线分别为患病率(x变量)为0.05、0.2、0.35、0.5、0.65时的曲线。
由图可知:当患病率小于0.35时,分组检测的方法更有效。患病率大于0.35,分组检测的方法基本不能减少工作量