3種方法對(duì)中肋骨條藻Logistic生長(zhǎng)模型擬合的比較研究

孫鵬 孫源

摘要：為探索高效準(zhǔn)確的Logistic生長(zhǎng)模型擬合方法，對(duì)中肋骨條藻在不同濃度二價(jià)鎘[Cd（Ⅱ）]脅迫下的Logistic生長(zhǎng)曲線，分別用Excel、Origin和Matlab 3種軟件進(jìn)行擬合。結(jié)果表明，擬合后的R2大小順序?yàn)镸atlab>Origin>Excel;采用Excel非線性擬合中的四點(diǎn)法擬合效果優(yōu)于三點(diǎn)法。提示Matlab（用Levenberg-Marquardt算法進(jìn)行非線性化多次迭代）是3種軟件中對(duì)Logistic生長(zhǎng)曲線擬合的最優(yōu)方法。

關(guān)鍵詞：海洋硅藻;中肋骨條藻;Logistic生長(zhǎng)曲線;相關(guān)系數(shù)R2;Matlab;模型擬合;最優(yōu)方法篩選

中圖分類號(hào)： Q141;S184 ?文獻(xiàn)標(biāo)志碼： A ?文章編號(hào)：1002-1302（2019）02-0229-04

Logistic曲線（或Logistic函數(shù)）模型由比利時(shí)數(shù)學(xué)家Verhulst首次提出[1]，最初主要是描述和研究人口增長(zhǎng)的一種S形曲線。其方程為y=a1+b·e-kx，隨后Logistic曲線在各領(lǐng)域得到了廣泛應(yīng)用，如Pearl等首次將Logistic曲線應(yīng)用于生物領(lǐng)域[2];Fisk在經(jīng)濟(jì)學(xué)領(lǐng)域引入Logistic曲線[3];Anderson將Logistic曲線應(yīng)用到醫(yī)療診斷中[4]。如今，Logistic曲線幾乎替代了一切的S形曲線，成為描述種群生態(tài)、動(dòng)植物及微藻生長(zhǎng)變化的一條最為普遍的曲線[5]。

Logistic曲線的擬合就是確定方程中的3個(gè)參數(shù)值a、b、k，經(jīng)典方法包括線性和非線性方法。線性方法是將Logistic曲線方程轉(zhuǎn)為線性方程，然后采用最小二乘法讓殘差平方和達(dá)到極小值以對(duì)參數(shù)進(jìn)行無(wú)偏估計(jì)。非線性方法可直接用相關(guān)算法得到Logistic模型參數(shù)的最小二乘無(wú)偏估計(jì)[6];現(xiàn)采用的算法有萊文貝格-馬夸特方法（Levenberg-Marquardt）方法、數(shù)值方法、三次樣條插值函數(shù)法等[7]。上述參數(shù)估計(jì)的過(guò)程較復(fù)雜，須借助各種軟件及編程語(yǔ)言來(lái)進(jìn)行擬合分析，如統(tǒng)計(jì)學(xué)軟件SPSS、SAS，程序設(shè)計(jì)語(yǔ)言R、Python。目前熱門的機(jī)器學(xué)習(xí)算法也涉足非線性擬合領(lǐng)域（如支持向量機(jī)回歸[8]），而掌握這些軟件和語(yǔ)言對(duì)生物從業(yè)者較為困難。

中肋骨條藻（Skeletonema costatum）是常見(jiàn)的海洋硅藻，分布較廣泛。近年來(lái)，研究者多對(duì)海洋微藻重金屬脅迫下的生長(zhǎng)過(guò)程用Logistic方程進(jìn)行描述，進(jìn)而反映出藻生物量與生長(zhǎng)時(shí)間、重金屬濃度的關(guān)系[9]。本研究以中肋骨條藻為模式生物，分別采用Excel、Origin、Matlab 3種常用的作圖和計(jì)算軟件，對(duì)其在二價(jià)鎘[Cd（Ⅱ）]脅迫下的生長(zhǎng)用Logistic模型進(jìn)行擬合，對(duì)擬合方法和效果進(jìn)行介紹和比較，為生命科學(xué)從業(yè)者進(jìn)行Logistic曲線的擬合、參數(shù)估計(jì)和數(shù)據(jù)處理提供借鑒和參考。

1 材料與方法

1.1 微藻培養(yǎng)及毒性試驗(yàn)

中肋骨條藻接種自中國(guó)海洋大學(xué)海洋生命學(xué)院，在曲阜師范大學(xué)生命科學(xué)學(xué)院環(huán)境生態(tài)學(xué)實(shí)驗(yàn)室采用Guillard F/2（f/2）培養(yǎng)基培養(yǎng)，培養(yǎng)條件為溫度（20±1） ℃，光源為白色日光燈，光照強(qiáng)度為3 000 lx，12 h光照12 h黑暗培養(yǎng);各步驟均進(jìn)行滅菌處理。

取20 mL藻液于250 mL f/2培養(yǎng)基中，微藻起始細(xì)胞密度均為104個(gè)/mL，加入4種不同濃度梯度（0、0.1、1.0、5.0 mg/L）的Cd（Ⅱ）試驗(yàn)，每個(gè)濃度下均設(shè)平行樣3個(gè)，采用與微藻培養(yǎng)相同的培養(yǎng)條件，每隔24 h取10 mL藻液，用Coulter計(jì)數(shù)儀測(cè)定藻液細(xì)胞密度，重復(fù)測(cè)定3次[10]，連續(xù)測(cè)定7 d。

1.2 數(shù)據(jù)處理軟件及系統(tǒng)

軟件：Microsoft Excel 2003、Origin 8.0、Matlab R2010b;系統(tǒng)：Windows7×86 with Service Pack1。

1.3 微藻生長(zhǎng)擬合

“1.1”節(jié)中的微藻生長(zhǎng)曲線用“1.2”節(jié)中的3種軟件進(jìn)行Logistic回歸擬合。

1.3.1 Excel對(duì)Logistic曲線的線性化擬合

利用公式（5）和a值可以先求出y′，利用Excel對(duì)y′和x進(jìn)行線性回歸分析，擬合后可得到b和k值。在Excel 2003中的操作如下：

（1）[工具]→[加載宏]→勾選“分析工具庫(kù)”。（2）[工具]→[數(shù)據(jù)分析]→[回歸]→彈出回歸對(duì)話框→單擊“Y值輸入?yún)^(qū)域（Y）”輸入框右邊的折疊按鈕，選取y′，單擊“X值輸入?yún)^(qū)域（X）”輸入框右邊的折疊按鈕，選取x的數(shù)據(jù);→在“輸出選項(xiàng)”框中選擇“輸出區(qū)域”單選框，點(diǎn)擊右邊的折疊按鈕，任選一輸出區(qū)域;→勾選“殘差”框中的“線性擬合”;→單擊[回歸]對(duì)話框中的[確定]，顯示回歸結(jié)果，并繪得散點(diǎn)圖。（3）利用散點(diǎn)圖，添加趨勢(shì)線可得到線性回歸公式。

1.3.2 Origin進(jìn)行Logistic曲線的非線性回歸擬合 （1）打開(kāi)Origin8.0，點(diǎn)選[column]→[add new columns]輸入須要增加的數(shù)據(jù)列數(shù)（本試驗(yàn)3列）。（2）按列輸入數(shù)據(jù)后，選中所有數(shù)據(jù)，菜單欄中選擇[Plot]→[Symbol]→[Scatter]繪制散點(diǎn)圖（圖1-a圖）。（3）關(guān)閉坐標(biāo)軸設(shè)置對(duì)話框，點(diǎn)選1-選擇黑色方框數(shù)據(jù)組b進(jìn)行擬合→菜單欄中選擇[Analysis]→[Fitting]→[Nonlinear Curve Fit]→[Open Dialogue]打開(kāi)對(duì)話框;按步驟點(diǎn)選2-Category-Growth/Sigmonial;3-Function-Slogistic3;4-Fit，即完成擬合和作圖（圖1-b），同時(shí)給出參數(shù)和統(tǒng)計(jì)處理數(shù)據(jù)。

1.3.3 用Matlab進(jìn)行Logistic曲線的非線性回歸擬合 打開(kāi)Matlab R2010b，用簡(jiǎn)單編程語(yǔ)言和自帶的控件工具箱對(duì)Logistic曲線快速擬合：（1）在Command Window對(duì)話框中輸入原始數(shù)據(jù)：橫坐標(biāo)x以及y1、y2、y3、y4這4條生長(zhǎng)曲線的數(shù)據(jù)。（2）依次點(diǎn)擊軟件左下角：[Start]→[Toolboxes]→[Curve Fitting]→[Curve Fitting Tool（cftool）]，彈出對(duì)話框“Curve Fitting Tool”。（3）在Curve Fitting Tool對(duì)話框中點(diǎn)擊Data→彈出Data對(duì)話框→選擇：X Data：x; Y Data：y1;→點(diǎn)擊Create data set按鈕可繪制散點(diǎn)圖;（4）在Curve Fitting Tool對(duì)話框中點(diǎn)擊Fitting，彈出Fitting對(duì)話框→在Fitting對(duì)話框中點(diǎn)選New fit，Type of fit：選擇Custom Equations，點(diǎn)擊New按

鈕→彈出New Custom Equation對(duì)話框，單擊選擇General Equations選項(xiàng)→在General Equations選項(xiàng)中輸入Equation：y=a/[1+b·exp（-k·x）]→點(diǎn)擊OK（圖2-a）。（5）回到Fitting選項(xiàng)卡，點(diǎn)擊Fit options→Algorithms：Levenberg-Marquardt算法→Close→Apply，即可得擬合后的Logistic曲線圖（圖2-b）。

2 結(jié)果與分析

2.1 3種方法繪制的不同濃度Cd（Ⅱ）下中肋骨條藻的S形生長(zhǎng)曲線

由圖3至圖5可得，（1）圖3中用Excel作出的曲線只能是連接各試驗(yàn)數(shù)據(jù)點(diǎn)的平滑曲線，無(wú)法擬合作圖;（2）圖5中Matlab所得曲線與各數(shù)據(jù)點(diǎn)的吻合程度明顯高于圖4，說(shuō)明Matlab對(duì)生長(zhǎng)曲線的擬合較優(yōu)。

2.2 3種擬合方法對(duì)微藻生長(zhǎng)曲線的擬合效果比較

由表1至表4可得，3種軟件擬合后的R2值整體大小順序?yàn)镸atlab>Origin>Excel，其中Excel方法中四點(diǎn)法R2值均大于三點(diǎn)法。數(shù)據(jù)同時(shí)顯示，四點(diǎn)法的參數(shù)值與非線性方法更為接近。在非線性擬合中，Matlab多次迭代的R2值均高于Origin的單次迭代。

3 結(jié)論與討論

Logistic曲線作為一種廣義的線性回歸模型，它的擬合方法可分為線性和非線性法。本研究3種軟件中，Excel對(duì)應(yīng)的是Logistic曲線的線性方法擬合;Origin和Matlab對(duì)應(yīng)的是曲線的非線性回歸擬合，其中Origin采用單次迭代，而Matlab是多次迭代。曲線擬合精度的高低都可用R2的大小來(lái)判定。由試驗(yàn)結(jié)果可知，3種軟件對(duì)Logistic的擬合精度均較高，R2值均接近于1[表4顯示Cd（Ⅱ）濃度較大時(shí)，對(duì)微藻的生長(zhǎng)會(huì)產(chǎn)生較大影響，R2值減小]，這表明3種軟件所獲得的Logistic回歸方程均與實(shí)測(cè)數(shù)據(jù)擬合很好。同時(shí)結(jié)果提示，（1）Excel可以對(duì)Logistic曲線模型的3個(gè)參數(shù)進(jìn)行估計(jì)，但不能作出擬合曲線圖，只能得到平滑曲線圖。（2）Excel進(jìn)行擬合時(shí)，四點(diǎn)法各參數(shù)估計(jì)量與非線性回歸最為接近，因此本研究中四點(diǎn)法比三點(diǎn)法擬合效果優(yōu)。這是由于本研究數(shù)據(jù)點(diǎn)數(shù)為8個(gè)，通常四點(diǎn)法較適合于試驗(yàn)數(shù)據(jù)序列數(shù)為偶數(shù)的情況，三點(diǎn)法適用于奇數(shù)的情況，這與以往報(bào)道[13]一致，但三點(diǎn)法具有簡(jiǎn)單快速的優(yōu)勢(shì)。（3）Excel擬合中所采用的線性擬合方法雖然簡(jiǎn)便，但對(duì)Logistic函數(shù)進(jìn)行線性變換后會(huì)對(duì)估計(jì)參數(shù)的性質(zhì)產(chǎn)生影響，如不再具有無(wú)偏性;所以其擬合精度不如Origin和Matlab所采用的非線性擬合。（4）在非線性擬合中，多次迭代會(huì)提高R2，因此Matlab的擬合效果優(yōu)于Origin。（5）綜合數(shù)據(jù)可得，不同擬合方法擬合后的R2大小順序?yàn)榉蔷€性>四點(diǎn)線性>三點(diǎn)線性，3種軟件擬合后的R2值為Matlab>Origin>Excel，其中，Matlab采用Levenberg-Marquardt算法進(jìn)行非線性化多次迭代擬合，提示Matlab是3種軟件中對(duì)Logistic生長(zhǎng)曲線擬合的最佳方法。

Logistic曲線擬合的優(yōu)度檢驗(yàn)除了R2，還有如卡方值、標(biāo)準(zhǔn)偏差等其他統(tǒng)計(jì)指標(biāo)[14]，除了本研究的3種軟件，對(duì)Logistic的擬合還可借助其他各種計(jì)算機(jī)語(yǔ)言。因此各種不同軟件、方法對(duì)Logistic擬合的比較還有更多的探討空間。

參考文獻(xiàn)：

[1]Verhulst P J. Notice sur la loi que la population suit dans sons accroissement[J]. Corr Math et Physique，1838，10：113-121.

[2]Pearl R，Reed L J. The Logistic Curve and the Census Count of 1930[J]. Science，1930，72（1868）：399-401.

[3]Fisk P R. The Graduation of income distributions[J]. Econometrica，1961，29（2）：171-185.

[4]Anderson J A. Sepatate sample Logistic discrimination[J]. Biometrika，1972，59：19-35.

[5]王福林，王吉權(quán). 生長(zhǎng)曲線參數(shù)估計(jì)的一種新方法——優(yōu)化回歸組合法[J]. 生物數(shù)學(xué)學(xué)報(bào)，2007，22（3）：533-538.

[6]陳 冬，程維虎. 利用樣本分位數(shù)的Logistic總體分布參數(shù)的近似最佳線性無(wú)偏估計(jì)[J]. 應(yīng)用數(shù)學(xué)學(xué)報(bào)，2005，28（2）：325-332.

[7]許小勇，鐘太勇. 三次樣條插值函數(shù)的構(gòu)造與Matlab實(shí)現(xiàn)[J]. 兵工自動(dòng)化，2006，25（11）：76-78.

[8]Zhang L J，Song Z J，Pan X D，et al. Comparison of v-support vector regression and Logistic equation for descriptive modeling of Lactobacillus plantarum growth[J]. African Journal of Biotechnology，2011，10（32）：6162-6171.

[9]欒紅艷，趙衛(wèi)紅，苗 輝. Cd2+脅迫下中肋骨條藻細(xì)胞內(nèi)多胺的生理響應(yīng)[J]. 中國(guó)環(huán)境科學(xué)，2015，35（5）：1487-1494.

[10]Torres E，Mera R，Herrero C，et al. Isotherm studies for the determination of Cd（Ⅱ） ions removal capacity in living biomass of a microalga with high tolerance to cadmium toxicity[J]. Environmental Science & Pollution Research，2014，21（22）：12616-12628.

[11]殷祚云. Logistic曲線擬合方法研究[J]. 數(shù)理統(tǒng)計(jì)與管理，2002，1（21）：41-46.

[12]胡文冉，范 玲，田曉莉，等. Excel在Logistic曲線擬合中的應(yīng)用[J]. 農(nóng)業(yè)網(wǎng)絡(luò)信息，2013（3）：14-16.

[13]范國(guó)兵. 一種估計(jì)Logistic模型參數(shù)的方法及應(yīng)用實(shí)例[J]. 經(jīng)濟(jì)數(shù)學(xué)，2010，27（1）：105-110.

[14]梁玉麟，趙建東. 基于logistic方程的供應(yīng)商-采購(gòu)商合作關(guān)系的經(jīng)濟(jì)模型研究[J]. 魯東大學(xué)學(xué)報(bào)（自然科學(xué)版），2015，31（3）：200-204，213.