楊益民，付必謙

（1.北京工商大學(xué) 理學(xué)院，北京，100084；2.首都師范大學(xué) 生命科學(xué)學(xué)院，北京100084）

0 引言

是描述有限環(huán)境條件下生物種群S-型增長的最基本和最常用的模型，在種群數(shù)量動態(tài)研究中占有十分重要的地位。

在利用式（1）定量研究種群動態(tài)時，必須依據(jù)實驗觀測數(shù)據(jù)估計式中的K、r和a。較早提出且目前仍然十分常用的方法是將式（1）移項取對數(shù)后，轉(zhuǎn)化為線性方程

先根據(jù)原始數(shù)據(jù)和經(jīng)驗預(yù)先給定K，然后再利用一元線性回歸估計r和a[1～3]。后來一些學(xué)者認(rèn)為線性化方法主觀色彩較濃或遺漏“最優(yōu)”參數(shù)估計值的可能性較大，提出采用一些非線性最優(yōu)化方法如麥夸法和單純形加速法等[4，5]，也有作者對利用數(shù)值方法擬合進(jìn)行了一些探討[6]。但是，無論是利用線性化還是非線性化方法（包括直接調(diào)用統(tǒng)計軟件包），在Logistic模型參數(shù)估計中均存在一些共同的需要特別重視的問題。此外，對現(xiàn)有線性化方法進(jìn)行改進(jìn)，探討一種既簡便易行又客觀準(zhǔn)確的線性化方法，對于缺乏較深數(shù)學(xué)基礎(chǔ)的實際工作者仍不失一種較好的選擇。

本文對目前Logistic模型參數(shù)估計中的問題進(jìn)行了分析，并對線性化方法進(jìn)行了再探討，提出了一種新的線性化方法——“最佳”觀測點數(shù)定向搜索法。將該方法應(yīng)用于實例的分析，獲得了滿意的擬合效果。

1 Logistic增長模型中a的估計方法

對于Logistic增長模型的擬合，絕大多數(shù)作者均將式（1）中的a、K和r作為3個獨立的參數(shù)進(jìn)行估計。在統(tǒng)計軟件包相關(guān)教程中，Verhulst模型的表達(dá)式也為y=b1/(1+b3exp(-b2x))，其中bi（i=1，2，3）為待估計的參數(shù)[7]。利用這些方法估計參數(shù)，計算得到的理論初值常與實測初值不等，甚至差異很大。雖然唐啟義等[8]曾對上述作法提出過質(zhì)疑，但似乎并未引起重視。

因此，Logistic模型參數(shù)估計所必須解決的一個重要問題是：能將a作為與K和r一樣的獨立的參數(shù)估計嗎？經(jīng)過擬合得到的理論初值是否應(yīng)與實測初值N0相同？

為探討上述問題，將Logistic增長模型：

則待擬合的方程應(yīng)為式（6）而非式（2）。

1 利用線性回歸分析估計參數(shù)K和r的新方法：“最佳”觀測點數(shù)定向搜索法

關(guān)于參數(shù)K和r的估計，許多線性化方法是首先估計一個K值，然后將其代入式（2）估計r，因此r的估計在很大程度上依賴于K的估計精度并受回歸分析采用的觀測點數(shù)目的影響。若K值估計精度不夠高或不確定性較大或選擇的觀測點數(shù)目不恰當(dāng)，均可能嚴(yán)重影響最終擬合結(jié)果。萬昌秀和梁中宇[3]提出的枚舉選優(yōu)法將種群達(dá)到平衡后的最大和最小觀測值作為K的上、下界，分別將該范圍內(nèi)的每個整數(shù)作為K值，并對每個K值均分別從全部觀測值中選取前m個觀測值（m=4，5，6，…，n；n為觀測值總數(shù)）進(jìn)行回歸分析估計可能的r值，然后逐一比較各組參數(shù)估計值對應(yīng)的SSe，以SSe最小為標(biāo)準(zhǔn)確定最終參數(shù)估計值。該方法考慮全面、穩(wěn)妥可靠，明顯提高了擬合優(yōu)度；但因計算量大（種群數(shù)量很大和觀測數(shù)據(jù)很多的實驗數(shù)據(jù)尤其如此）而需編制專門的計算機程序，并難以處理非整數(shù)型數(shù)據(jù)或種群未達(dá)平衡狀態(tài)便終止實驗的不完整的觀測數(shù)據(jù)。

針對上述問題，在提出利用式（6）作為擬合方程的基礎(chǔ)上，本文提出一種新的線性化方法——“最佳”觀測點數(shù)定向搜索法。該方法的要點為：（1）首先確定適當(dāng)?shù)腒值初始搜索區(qū)間；（2）利用該區(qū)間端點和中值確定一個在整個K值搜索過程中進(jìn)行線性回歸分析時均可采用的適宜的觀測點數(shù)目，從而不必對每個K值均選取前m個觀測值進(jìn)行分析；（3）利用黃金分割法定向搜索K的最優(yōu)估計值。

1.1 K值初始搜索區(qū)間[Ka，Kb]的確定

1.1.1 完整實驗數(shù)據(jù)初始搜索區(qū)間的確定

對于實驗數(shù)據(jù)相對完整、有多個觀測值圍繞平衡數(shù)量波動的情況，一般來講，參數(shù)K的最優(yōu)估計值K*應(yīng)在這些觀測值的平均值（Km）附近，因此可直接將Km作為初始搜索區(qū)間的中值Kmid來確定初始搜索區(qū)間[Ka，Kb]（Kmid=(Ka+Kb)/2）。若在這些數(shù)據(jù)中存在個別觀測值明顯偏大或偏小，可去除這些觀測值后再計算Km。初始搜索區(qū)間應(yīng)稍大一些，最好包括圍繞平衡數(shù)量波動的所有觀測值（異常值除外）。

1.1.2 不完整實驗數(shù)據(jù)初始搜索區(qū)間的確定

若種群未達(dá)平衡狀態(tài)便終止了實驗，則無法依據(jù)上述方法確定Kmid和初始搜索區(qū)間，而需預(yù)先估計一個搜索區(qū)間，然后再行調(diào)整。

（1）初始搜索區(qū)間的預(yù)估。對于增長曲線上端已經(jīng)明顯變緩的實驗數(shù)據(jù)，可直接通過目測預(yù)估初始搜索區(qū)間[K'a，K'b]，計算中值K'mid。若實驗數(shù)據(jù)明顯不完整，難以估計上漸近線的大致位置，則可隨機選取3組等時間間距觀測數(shù)據(jù)依據(jù)Pearl和Reed的三點法十分方便地估計3個K值，計算其均值Km，以Km作為區(qū)間中值K'mid初步估計初始搜索區(qū)間[K'a，K'b]。初始搜索區(qū)間應(yīng)包括上述3個K值并向上、下有一定延伸；對于增長曲線上端已出現(xiàn)轉(zhuǎn)折的實驗數(shù)據(jù)，還應(yīng)包括轉(zhuǎn)折段的觀測值，以減小遺漏K*的可能性。

1.2 回歸分析采用的“最佳”觀測點數(shù)（M★）的確定

在K值搜索過程中，若初始搜索區(qū)間[Ka，Kb]設(shè)置比較合理，K*一般應(yīng)位于區(qū)間中值Kmid附近，搜索過程中（尤其后期）大多數(shù)探測點的K值往往也距Kmid相對較近，因此可利用區(qū)間端點和中值探測一個在整個K值搜索過程中進(jìn)行回歸分析時均可采用的適宜觀測點數(shù)（“最佳”觀測點數(shù)M*）。具體方法如下：

（1）對于初始搜索區(qū)間[Ka，Kb]，分別根據(jù)式（6）計算端點Ka、Kb及中值Kmid對應(yīng)的Y值，在同一圖中作3個K值對應(yīng)的Y-t折線，然后在圖中找到符合下列條件的一點：在該點和該點之前，不僅3條Y-t折線的線性關(guān)系均較強且比較類似，而且包括的觀測點數(shù)目盡可能多，使依據(jù)這些觀測數(shù)據(jù)估計r時能夠利用更多實驗數(shù)據(jù)的信息和獲得更高的擬合優(yōu)度?？梢钥吹?，該點一般對應(yīng)于3條Y-t折線呈明顯“Y”形分離的分叉點（將其對應(yīng)的觀測點數(shù)記為ms），但有時也可能為分叉點之前或之后的一點，因此比較適宜的觀測點數(shù)分別為ms-1、ms和ms+1。

本方法的關(guān)鍵在于以Ka、Kb及Kmid對應(yīng)的3條Y-t折線確定ms，因此可稱為“三線定點法”。

“三線定點法”的合理性證明：

在初始搜索區(qū)間[Ka，Kb]中，設(shè) Kc=Ka+ΔK(0≤ΔK≤Kb-Ka)，則對于任一Kc，均有

因此，在Y-t圖中，隨著觀測點數(shù)目的變化（反映N的變化），不僅特定K值對應(yīng)的Yc值的變化趨勢、而且不同K值對應(yīng)的Yc值彼此之間的大小關(guān)系均可能發(fā)生明顯改變，由此必然影響參數(shù)r的估計結(jié)果。隨著N的增加，各個Kc=Ka+ΔK對應(yīng)的Y-t折線將按照ΔK的大小依次向下發(fā)生轉(zhuǎn)折，直至Kc大于全部N值。在區(qū)間[Ka，Kb]中，Ka最小，因此在N＜Ka的數(shù)據(jù)范圍內(nèi)，各個K值對應(yīng)的Y-t折線不僅均處于上升段而且位于兩端點對應(yīng)的Y-t折線之間。若以區(qū)間中值Kmid對應(yīng)的Y-t折線作為參照線，則在Ka、Kb和Kmid對應(yīng)的3條Y-t折線均呈較強線性關(guān)系的觀測點范圍內(nèi)，其他任一Kc值對應(yīng)的Y-t折線也應(yīng)保持類似關(guān)系，因此可在整個K值搜索過程中以同樣的觀測點數(shù)目進(jìn)行回歸分析估計r。

以Gause（1934）雙小核草履蟲種群增長實驗（添加一環(huán)細(xì)菌）[9]為例。該實驗數(shù)據(jù)完整，有多個觀測值圍繞平衡數(shù)量上下波動，因此可直接以這些觀測值的平均值（Km＝448）作為中值Kmid確定初始搜索區(qū)間「418，478」。為了更加清晰地反映不同K值對應(yīng)的Y-t折線之間的關(guān)系，避免實測值在平衡數(shù)量附近的波動對圖形的影響，首先以唐啟義等[16]的預(yù)測值代替實測值進(jìn)行分析（圖1左圖）。在預(yù)測值中，N7=409＜Ka，N8=429＞Ka，Nmax=N16=439。由圖可見，在前8個觀測點范圍內(nèi)，所有Kc=Ka+ΔK(0＜ΔK≤Kb-Ka)對應(yīng)的Yc-t折線均位于Ya-t折線之下且與Ya的偏離程度隨△K的增大而增大。在該范圍之后，對于所有KC＜Nmax=439，其對應(yīng)的Yc-t折線均隨N的增加依次向下轉(zhuǎn)折，最終的Yc值隨ΔK的增大而增大。而對于其他Kc值，因KC＞Nmax，其對應(yīng)的Yc值單調(diào)增加直至穩(wěn)定，且隨ΔK的增大而減小。各條Y-t折線明顯“Y”形分叉點對應(yīng)的觀測點數(shù)ms=7。

若直接以實驗觀測值計算Yc，則得圖1右圖。實驗觀測值N8=396＜Ka，N9=443＞Ka。由圖可見，盡管該圖右部各Y-t折線波動劇烈，但在前9個觀測點范圍內(nèi)，Y-t折線的排列卻十分有序，Yc與Ya的偏離程度隨ΔK的增大而增大。各條Y-t折線明顯“Y”形分叉點對應(yīng)的觀測點數(shù)也為ms=7。

圖1 “三線定點法”示意圖。依據(jù)Gause（1934）雙小核草履蟲種群增長數(shù)據(jù)（添加一環(huán)細(xì)菌）[9]作圖。左：利用唐啟義等預(yù)測值[8]計算Y；右：利用Gause實測值計算Y。

1.3 K的“最優(yōu)”估計值（K★）的定向搜索

黃金分割法是運籌學(xué)中單變量非線性最優(yōu)化的重要方法[18]，該法具有搜索方向明確、收斂速度快、計算簡單且計算量小的優(yōu)點，而且適用于整數(shù)型或非整數(shù)型觀測數(shù)據(jù)。利用該法搜索K*的具體作法如下：

（1）對于初 始搜索區(qū)間[Ka，Kb]，搜索區(qū)間 長度L=Kb-Ka。在搜索區(qū)間內(nèi)距兩個端點各0.382 L處對稱插入兩個內(nèi)插點Kc和Kd，下內(nèi)插點Kc=Ka+0.382L，上內(nèi)插點Kd=Kb-0.382L=Ka+0.618L。分別將Kc和Kd代入式（6），選取前M*個觀測數(shù)據(jù)進(jìn)行回歸分析估計r，然后計算對應(yīng)的SSe。

（2）比較兩內(nèi)插點對應(yīng)的SSe，若下內(nèi)插點Kc的SSe較大，則舍棄下端點Ka，將Kc作為新區(qū)間的端點，并在上內(nèi)插點Kd與上端點Kb之間插入一個新的上內(nèi)插點；反之，若上內(nèi)插點Kd的SSe較大，則舍棄上端點Kb，將Kd作為新區(qū)間端點，并在Kc與Ka之間插入新的下內(nèi)插點。仍按照（?。┲械姆椒ㄓ嬎阈聝?nèi)插點的K值（此時L應(yīng)為縮小后的搜索區(qū)間長度）。

（3）多次重復(fù)步驟（ⅱ），通過不斷舍棄端點來逐漸縮小搜索范圍（第n次搜索后，搜索長度將縮小為初始長度L的0.618n倍），同時不斷插入新內(nèi)插點以加大對重點區(qū)域的搜索密度。當(dāng)SSe下降十分緩慢或達(dá)到實驗要求的精度時，即可將此時的K值作為“最優(yōu)”估計值K*。

1.4 方法小結(jié)

綜上所述，“最佳”觀測點數(shù)定向搜索法的主要分析步驟如下：

（1）完整實驗數(shù)據(jù)的擬合。

①實驗觀測數(shù)據(jù)的預(yù)處理。若原始實驗數(shù)據(jù)的初始時間t0≠0，則以t'=t-t0代替t。

②選取圍繞平衡數(shù)量波動的觀測值計算其平均值Km，以Km作為中值Kmid，確定初始搜索區(qū)間[Ka，Kb]。。

③利用“三線定點法”確定一個在整個K值搜索過程中進(jìn)行線性回歸時均可采用的“最佳”觀測點數(shù)M*。

④利用黃金分割法定向搜索參數(shù)K的“最優(yōu)”估計值K*和相應(yīng)的r值。

⑤此外，搜索到K*值之后，還可對其進(jìn)一步調(diào)整回歸分析采用的觀測點數(shù)目，以提高參數(shù)r的估計精度、獲得更好的擬合效果。

⑥將K*和最終獲得的r估計值r*代入式（6）或式（1）（積分常數(shù)a據(jù)K*和實測初值N0由式（5）解出），即可獲得Logistic增長模型。

（2）不完整實驗數(shù)據(jù)的擬合。

若實驗數(shù)據(jù)不完整，則需首先預(yù)估一個K值搜索區(qū)間（具體方法見1.1.2步驟（ⅰ）），作該區(qū)間端點和中值對應(yīng)的Y-t圖，確定ms；然后分別將端點和中值代入式（6）估計r（因兩個端點的K值往往相差較大，一般直接由前ms個觀測數(shù)據(jù)進(jìn)行回歸分析即可），計算3個K值對應(yīng)的SSe。若一側(cè)端點的SSe最低，則以該端點作為新區(qū)間的中值，將搜索區(qū)間向該側(cè)平移。重復(fù)上述調(diào)整過程，直至中值的SSe最低，最終確定初始搜索區(qū)間[Ka，Kb]。

其他步驟與上述完整實驗數(shù)據(jù)的分析相同。

2 實例分析及結(jié)果

利用“最佳”觀測點數(shù)定向搜索法對7個實例進(jìn)行了分析。其中，Gause（1934）草履蟲種群增長實驗[9]是種群生態(tài)學(xué)中最為經(jīng)典和重要的實驗。

2.1 完整實驗數(shù)據(jù)的分析過程

仍以Gause（1934）雙小核草履蟲種群增長實驗（添加一環(huán)細(xì)菌）[9]為例。由前面分析可知，本例初始搜索區(qū)間[418，478]，Kmid＝448，ms=7。取Kmid＝448，分別以前6、7、8個觀測數(shù)據(jù)估計r。3個r值對應(yīng)的SSe分別為17436、14504和14970。ms=7對應(yīng)的SSe最小，因此確定M*=7。以前7個觀測數(shù)據(jù)進(jìn)行回歸分析估計r，經(jīng)8次搜索（表1），得K*=438，r*=1.1532。

表1 Gause（1934）雙小核草履蟲實驗數(shù)據(jù)（添加一環(huán)細(xì)菌）[9]K值搜索過程

2.2 不完整實驗數(shù)據(jù)初始搜索區(qū)間的確定

以王振中和林孔勛花生銹病顯癥動態(tài)試驗（Ⅰ）前19天的實驗數(shù)據(jù)[12]為例。該例觀測時間為第9-19天，相應(yīng)的觀測值（孢子堆數(shù)目）依次為79、126、257、625、1856、4000、7389、12133、16057、18898和22361。首先將原觀測時間轉(zhuǎn)換為 t=0～10 ；然后分別選取t=（0，5，10）、（1，5，9）和（2，5，8）對應(yīng)的觀測值，按照三點法[1]計算K值，得K=21564、17843和19204。將3個K值的均值Km=19537作為K'mid，分別向上、下各延伸4000預(yù)估初始搜索區(qū)間，得[15537，23537]。分別作K=15537、19537和23537對應(yīng)的Y-t圖，3條折線的明顯“Y”形分叉點對應(yīng)于ms=7。選取前7個觀測數(shù)據(jù)進(jìn)行回歸分析估計3個K值對應(yīng)的r值，計算SSe，分別為107790143、33024912和4105099。由于上端點K=23537對應(yīng)的SSe最低，因此將其作為新區(qū)間的中值，將初始搜索區(qū)間平移至[19537，27537]。對新區(qū)間進(jìn)行類似分析，得ms=8；選取前8個觀測數(shù)據(jù)估計r，以中值對應(yīng)的SSe最小，因此確定初始搜索區(qū)間為[19537，27537]。

2.3 個實例的分析結(jié)果

表2給出了7個實例的分析結(jié)果，相應(yīng)的Y-t折線圖分別見圖1-圖7。獲得“最優(yōu)”估計值K*后，均對其進(jìn)一步調(diào)整m值，以期提高參數(shù)r的估計精度，但最終均以m=M*對應(yīng)的 SSe最低。

表2 “最佳”觀測點數(shù)定向搜索法應(yīng)用實例分析

周賽花[13]誤將鄔祥光[4]第17天的實測值23.50看作28.50，因此擬合結(jié)果存在一定偏差。但為與周賽花提出的逐次加密

搜索法進(jìn)行比較，本文仍按其給出的數(shù)據(jù)進(jìn)行擬合。

圖2 Gause（1934）雙小核草履蟲（添加半環(huán)細(xì)菌）[17]Y-t圖

圖3 Gause（1934）大草履蟲（添加一環(huán)細(xì)菌）[17]Y-t圖

圖4 Gause（1934）大草履蟲（添加半環(huán)細(xì)菌）[17]Y-t圖

圖5 Carlson（1913）酵母種群[19]Y-t圖

圖6 稻癭蚊黃柄黑蜂[8]Y-t圖

圖7 花生銹病顯癥動態(tài)試驗Ⅰ（前19天數(shù)據(jù)）[7]Y-t圖

3 與其他作者擬合結(jié)果的比較

表3 Gause（1934）雙小核草履蟲種群增長實驗數(shù)據(jù)[17]擬合結(jié)果的比較＊

表4 Gause（1934）大草履蟲種群增長實驗數(shù)據(jù)[17]擬合結(jié)果的比較＊

表5 Carlson（1913）酵母種群增長實驗數(shù)據(jù)[19]擬合結(jié)果的比較＊

表6 稻癭蚊黃柄黑蜂種群增長數(shù)據(jù)[8]擬合結(jié)果的比較＊

表8 花生銹病顯癥動態(tài)試驗（Ⅰ）前19天實驗數(shù)據(jù)[7]擬合結(jié)果的比較＊

由上可見，與其他擬合方法比較，本文提出的“最佳”觀測點數(shù)定向搜索法具有如下優(yōu)點：

（1）不僅可以滿足Logistic增長模型的初值條件，保證理論初值與實測初值相等，使參數(shù)K和r的估計更為客觀準(zhǔn)確，而且能夠獲得理想的擬合效果，擬合優(yōu)度與唐啟義等[16]以及SPSS非線性回歸模塊（依據(jù)式（4））采用麥夸法的結(jié)果十分接近。

（2）每一步的分析結(jié)果都比較明確，最終擬合結(jié)果的不確定性小。

（3）可用于完整或不完整、整數(shù)型或非整數(shù)型實驗數(shù)據(jù)的擬合，適用范圍較之前的很多線性化方法明顯拓廣。

（4）原理簡單，計算量小，僅靠Excel即可輕松完成全部計算，很適于不熟悉統(tǒng)計軟件包的研究者利用。

[1]Andrewartha H G,Birch L G.The Distribution and Abundance of Animals[M].Chicago:University of Chicago Press,1954.

[2]考克斯.普通生態(tài)學(xué)實驗手冊[M].蔣有緒譯.北京：科學(xué)出版社，1979.

[3]萬昌秀，梁中宇.邏輯斯諦曲線的一種擬合方法[J].生態(tài)學(xué)報，1983，3（3）.

[4]王莽莽，李典謨.用麥夸方法最優(yōu)擬合邏輯斯諦曲線[J].生態(tài)學(xué)報，1986，6（2）.

[5]馬占山.單純形加速法擬合生態(tài)學(xué)中的非線性模型[J].生物數(shù)學(xué)學(xué)報，1992，7（2）.

[6]吳新元.邏輯斯諦曲線擬合的一種數(shù)值方法[J].生物數(shù)學(xué)學(xué)報，1990，5（1）.

[7]陸紋岱.SPSS for Windows統(tǒng)計分析[M].(第3版).北京：電子工業(yè)出版社，2006.

[8]唐啟義，胡國文，馮明光，胡陽.Logistic方程參數(shù)估計中的錯誤與修正[J].生物數(shù)學(xué)學(xué)報，1996，11（4）.

[9]Gause G F.The struggle for existence[M].Baltimore:Williams and Wilkins，1934.

[10]瓦格納.運籌學(xué)原理與應(yīng)用(第2版)[M].鄧三瑞等譯.北京：國防工業(yè)出版社，1992.

[11]Pearl R.The Growth of Populations[J].Quarterly Review of Biology.1927.

[12]王振中，林孔勛.邏輯斯諦曲線K值的四點式平均值估計法[J].生態(tài)學(xué)報，1987，7（3）.

[13]周賽花.邏輯斯諦方程中參數(shù)的估計[J].數(shù)理統(tǒng)計與管理，1992，11（5）.

[14]鄔祥光.昆蟲生態(tài)學(xué)的常用數(shù)學(xué)分析方法(修訂版)[M].北京：農(nóng)業(yè)出版社，1985.