標(biāo)記重捕法的模擬

標(biāo)記重捕法的模擬北京大學(xué)附屬中學(xué)肖然1.研究題目標(biāo)記重捕法，在被調(diào)查種群的生存環(huán)境中，捕獲一部分個(gè)體，將這些個(gè)體進(jìn)行標(biāo)志后再放回原來(lái)的環(huán)境，經(jīng)過(guò)一段時(shí)間后進(jìn)行重捕，根據(jù)重捕中標(biāo)志個(gè)體占總捕獲數(shù)的比例來(lái)估計(jì)該種群的數(shù)量。理論計(jì)算公式為：N=M×n/m。其中，N-估計(jì)種群數(shù)量，M-標(biāo)記數(shù)量，n-重捕數(shù)量，m-重捕個(gè)體中被標(biāo)記個(gè)體的數(shù)量。本文利用python編程工具，對(duì)標(biāo)記重捕法進(jìn)行了模擬仿真，驗(yàn)證了其正確性并進(jìn)行了簡(jiǎn)單的誤差分析。2.模擬分析為了科學(xué)地對(duì)標(biāo)記重捕法進(jìn)行模擬，我們將標(biāo)記重捕的整個(gè)過(guò)程分解為以下幾個(gè)過(guò)程：隨機(jī)產(chǎn)生所有種群個(gè)體的位置；合理地捕獲一些個(gè)體，并對(duì)捕獲的個(gè)體進(jìn)行標(biāo)記；重新捕獲一些個(gè)體；統(tǒng)計(jì)計(jì)算，根據(jù)標(biāo)記重捕原理N=M×n/m估算種群總數(shù)量，并與準(zhǔn)確值進(jìn)行對(duì)比。3.研究過(guò)程3.1隨機(jī)產(chǎn)生個(gè)體位置我們假設(shè)種群個(gè)體均勻分布在某個(gè)正方形的區(qū)域內(nèi)，不妨認(rèn)為此正方形的邊長(zhǎng)=2（否則可以根據(jù)真實(shí)值進(jìn)行縮放）。令(x,y)表示個(gè)體位置，利用numpy庫(kù)中的random.uniform(L,R)方法，將x,y設(shè)置為在[-1,1]范圍內(nèi)均勻分布的隨機(jī)數(shù)：importnumpyx=numpy.random.uniform(-1,1)#x,y設(shè)置為[-1,+1]范圍內(nèi)的隨機(jī)數(shù)y=numpy.random.uniform(-1,1)利用一個(gè)for循環(huán)，不斷地產(chǎn)生總共N只個(gè)體，并將它們的(x,y)坐標(biāo)記錄在兩個(gè)列表中。并利用matplotlib庫(kù)能夠?qū)㈦S機(jī)產(chǎn)生的(x,y)點(diǎn)畫在二維坐標(biāo)系中：importmatplotlib.pyplotaspltimportnumpyN=100xList=[]#定義2個(gè)空列表yList=[]foriinrange(N):x=numpy.random.uniform(-1,1)y=numpy.random.uniform(-1,1)xList.append(x)#將剛產(chǎn)生的x加入xListyList.append(y)#將剛產(chǎn)生的y加入yListplt.plot(xList,yList,'ro')plt.title('N=100')圖1.圖1.隨機(jī)產(chǎn)生的N只個(gè)體坐標(biāo)3.2捕獲一些個(gè)體并進(jìn)行標(biāo)記做一個(gè)簡(jiǎn)單抽象，將“捕獲”模擬為以一個(gè)特定概率（例如P=0.2）被選中，然后分兩組列表{xList1,yList1}/{xList2,yList2}，分別存放被/沒被捕獲的個(gè)體坐標(biāo)。模擬以特定概率發(fā)生某事，可以用uniform函數(shù)產(chǎn)生1個(gè)[0,1]內(nèi)隨機(jī)數(shù)，然后判斷其是否<P=0.2即可：foriinrange(N):x=numpy.random.uniform(-1,1)y=numpy.random.uniform(-1,1)isCap=numpy.random.uniform(0,1)<=0.2ifisCap:xList1.append(x)yList1.append(y)else:xList2.append(x)yList2.append(y)plt.plot(xList1,yList1,'ro')#標(biāo)記的為紅色plt.plot(xList2,yList2,'bo')#未標(biāo)記的為藍(lán)色圖2.圖2.捕獲一些個(gè)體進(jìn)行標(biāo)記plt.show()此時(shí)可以根據(jù)存放被捕獲個(gè)體坐標(biāo)的xList的大小計(jì)算出標(biāo)記數(shù)量M=len(xList1)。3.3重新捕獲一些個(gè)體“重捕”也可以模擬為以一個(gè)特定概率（例如P=0.2）被選中，這時(shí)種群個(gè)體會(huì)分為3類：重捕中捕獲的個(gè)體中被/沒被標(biāo)記過(guò)的，沒在重捕中捕獲的個(gè)體，我們分別用3組列表{xListA,yListA}/{xListB,yListB}/{xListC,yListC}存放以上3類個(gè)體的坐標(biāo)。重捕中捕獲的個(gè)體數(shù)量n=len(xListA)+len(xListB)，其中被標(biāo)記過(guò)的個(gè)體數(shù)量m=len(xListA)。圖圖3.重新捕獲一些個(gè)體，統(tǒng)計(jì)計(jì)算3.4統(tǒng)計(jì)計(jì)算設(shè)置真實(shí)種群數(shù)量N0=100，執(zhí)行python代碼，模擬結(jié)果為：M=20，n=29,m=7；根據(jù)理論計(jì)算公式N=M×n/m，計(jì)算出種群數(shù)量的估計(jì)值N=82.86，相比較真實(shí)值誤差較大；設(shè)置N0=10000，計(jì)算出種群數(shù)量的估計(jì)值N=9993.46，能夠發(fā)現(xiàn)相對(duì)誤差有所減小。嘗試不同的N0值，得到以下結(jié)果：真實(shí)種群數(shù)量N010010001000050000100000種群數(shù)量估計(jì)值N82.86918.429993.4649820.29100084.04相對(duì)誤差17.14%8.16%0.07%0.37%0.08%圖5.多次模擬取平均值圖圖5.多次模擬取平均值圖4.相對(duì)誤差跟隨N0的變化趨勢(shì)通過(guò)上述模擬分析，能夠發(fā)現(xiàn)選取N0=10000時(shí)，相對(duì)誤差已經(jīng)<0.1%，接著我們使用多次模擬取平均值的方法，進(jìn)一步減小誤差。importmatplotlib.pyplotaspltimportnumpyN=10000NMC=100#進(jìn)行100次模擬，畫出平均值的變化趨勢(shì)ans=[]ansmean=[]forkinrange(NMC):M=0m=0n=0foriinrange(N):ifnumpy.random.uniform(0,1)<0.2:M+=1foriinrange(M):ifnumpy.random.uniform(0,1)<0.3:m+=1n+=1foriinrange(N-M):ifnumpy.random.uniform(0,1)<0.3:n+=1print(k,M,n,m,M*n/m)ans.append(M*n/m)ansmean.append(sum(ans)/(k+1))plt.plot(ansmean)plt.show()4.研究結(jié)論通過(guò)對(duì)標(biāo)記重捕法在python中的模擬，驗(yàn)證了其正確性，并發(fā)現(xiàn)在種群數(shù)量較大時(shí)才會(huì)有可以接受的相對(duì)誤差。在這個(gè)過(guò)程中，利用了1）numpy中的random.uniform函數(shù)產(chǎn)生均勻分布的隨機(jī)數(shù)；2）列表List存放種群個(gè)體的坐標(biāo)；3）matplotlib中的plot函數(shù)進(jìn)行數(shù)據(jù)可視化。5.參考文獻(xiàn)[