劉 薇，常振海，張德生

(1.天水師范學(xué)院 數(shù)學(xué)與統(tǒng)計學(xué)院，甘肅 天水 741001；2.西安理工大學(xué) 理學(xué)院,西安 710054)

0 引言

糧食產(chǎn)量和糧食安全是人們一直關(guān)注的熱點(diǎn),學(xué)者們從不同角度討論了糧食問題。因為能夠用于統(tǒng)計分析的糧食年度數(shù)據(jù)較少,導(dǎo)致建立起的模型統(tǒng)計性質(zhì)缺乏可信度,bootstrap方法給我們提供了增加樣本量的便捷途徑,關(guān)于bootstrap方法的應(yīng)用研究近年來也熱度不減[1～4]。將bootstrap方法用于研究糧食產(chǎn)量的文獻(xiàn)所見不多,考慮到糧食產(chǎn)量前后幾十年增加較快,可能存在距離均值較遠(yuǎn)的數(shù)據(jù),本文結(jié)合較為穩(wěn)健的M-估計[5],建立起關(guān)于我國糧食產(chǎn)量的bootstrap法回歸模型。

1 基本知識與方法

考慮最簡單的回歸模型

在c=1.345時，它對正態(tài)數(shù)據(jù)的效率為95%。

(2)計算bootstrap法估計T*=β(S*)。

(3)重復(fù)(1)和(2)R次,本文取R=1999,這樣得到的R個T*r(r=1，…，R)。

(4)利用R個T*r得到T的近似分布,從而可以得到待估計參數(shù)β的點(diǎn)估計,區(qū)間估計等統(tǒng)計量。

很顯然,這種方法省去了對分布的假定及基于假定分布的統(tǒng)計推斷,利用重復(fù)抽樣來擴(kuò)大樣本量,建立起基于樣本的統(tǒng)計分布。

2 實證分析

樣本范圍是1980～2012年共33個樣本,初步選用變量為：

x1：農(nóng)業(yè)機(jī)械總動力(萬千瓦),

x2：受災(zāi)面積(千公頃),

x3：國家財政用于農(nóng)業(yè)的支出(億元),

x4：有效灌溉面積(千公頃),

x5：農(nóng)用化肥施用折純量(萬噸),

x6：農(nóng)村用電量(億千瓦小時),

x7：糧食作物播種面積(千公頃),

Y：糧食產(chǎn)量(萬噸),

數(shù)據(jù)來源于國家統(tǒng)計局網(wǎng)站。計算主要基于R語言。

自變量個數(shù)相對于樣本個數(shù)來說較多,都選入模型顯然不合適,為了選取對糧食產(chǎn)量有重要影響的因素,我們首先對這些自變量進(jìn)行了降維分析,結(jié)果見圖1。

圖1 因子分析碎石圖

圖2 最大方差法因子旋轉(zhuǎn)圖

從圖1能看出,提取兩個因子,方差累計貢獻(xiàn)率就達(dá)到了88.95%,這樣,在這7個自變量中抽取出兩個公因子,僅損失約11%的信息,效果還是很不錯的。為了能更加清晰地顯示提取的兩個因子及其實際代表意義,我們做了最大方差法因子旋轉(zhuǎn),結(jié)果見圖2。

從圖2中能看出,變量x1,x3,x4,x5,x6代表了第一因子,x2,x7代表了第二因子。我們結(jié)合這些變量的實際意義,很容易看出第一因子與糧食單位面積產(chǎn)量有關(guān),第二因子于糧食作物播種面積有關(guān)。進(jìn)一步分析,在第一因子中的變量,對糧食產(chǎn)量影響最大的因素是x5,在第二因子的變量中則是x7。因此,我們分別以這兩個變量為自變量,建立回歸模型,同時為了便于比較自變量對因變量的貢獻(xiàn)大小,這里首先對數(shù)據(jù)進(jìn)行了中心化處理,處理后的變量分別稱為zx5,zx7,zy,結(jié)果見表1(前4列)。

表1 M-估計結(jié)果

從表1中的t值能看出,系數(shù)是顯著的,不過這個推斷是基于大樣本去推斷的,很顯然,穩(wěn)妥起見,需要擴(kuò)大樣本量,采用bootstrap方法,結(jié)果見表1最后三列,這里的偏差和標(biāo)準(zhǔn)誤差分別是

和表1中M-估計方法相比,bootstrap M-估計方法的系數(shù)幾乎沒什么改變,偏差很小,標(biāo)準(zhǔn)誤略有增大,這個可能和離群值有關(guān)系,后面再討論離群值的問題。這些bootstrap方法的估計系數(shù)分布見圖3。

圖3 bootstrap法的M-估計系數(shù)分布

從圖3能看出,zx5系數(shù)接近于正態(tài)分布,尾部稍有點(diǎn)厚,zx7系數(shù)明顯右偏,和正態(tài)分布相去甚遠(yuǎn)。在預(yù)測中這兩個系數(shù)的關(guān)系可用圖4來顯示。

圖4 兩個系數(shù)的聯(lián)合分布散點(diǎn)圖

從圖4能看出,這兩個系數(shù)近似于不相關(guān),說明在對zy預(yù)測時它們相對獨(dú)立,這和開始我們采用降維的方法選取的變量結(jié)論是一致的。有了它們各自的分布和聯(lián)合分布,我們就可以計算這些系數(shù)的M-估計區(qū)間,結(jié)果見表2。

表2 bootstrap法的M-估計區(qū)間

從表2能看出,系數(shù)zx5的三種區(qū)間結(jié)果相差不是很大,和該系數(shù)的分布近似于正態(tài)分布有關(guān)系,但系數(shù)zx7的則相差較大,圖3中顯示了其分布左尾較長,這個結(jié)果和圖3中呈現(xiàn)的分布相一致,比較來看,顯然BCa區(qū)間的結(jié)果更加可信。

這樣,建立起的基于bootstrap法和M-估計方法的回歸模型為

比較起見,我們將基于普通最小二乘方法的回歸結(jié)果列出,見式(5)。

比較式(4)和(5),M-估計方法調(diào)低了農(nóng)用化肥施用折純量對糧食產(chǎn)量的影響,同時增加了糧食作物播種面積的影響。式(4)中系數(shù)的t值更大,標(biāo)準(zhǔn)誤更小,說明M-估計更可靠。從回歸結(jié)果看,在固定一個自變量的情形下,農(nóng)用化肥施用折純量和糧食作物播種面積每增加一個單位,糧食產(chǎn)量將分別增加1.0645和0.3612個單位,說明化肥在糧食種植中的作用較大,在我國經(jīng)濟(jì)快速發(fā)展的形勢下,播種面積增加較為困難,重視化肥的使用將有助于糧食產(chǎn)量的增加。模型擬合的殘差見圖5。

圖5 回歸擬合的殘差

從圖5能看出,1980年、1981年和2003年的數(shù)據(jù)明顯偏離樣本均值,說明原始數(shù)據(jù)中可能存在異常值,為了直觀顯示每個觀測點(diǎn)對系數(shù)的影響,采用Jackknife-after-Bootstrap方法,結(jié)果見（圖6）。

圖6 Jackknife-after-Bootstrap方法散點(diǎn)圖

從圖6能看出,觀測點(diǎn)1(1980年),2(1981年)增大了系數(shù)zx5,減小了系數(shù)zx7,而點(diǎn)33(2012年)則使兩個系數(shù)均增大。對照原始數(shù)據(jù)去看,因為我國糧食產(chǎn)量在近30年來增加較快,所以從數(shù)據(jù)上就使得1980年和1980年的糧食產(chǎn)量數(shù)據(jù)顯得異常,而第24個數(shù)據(jù)(2003年)是近些年來糧食產(chǎn)量唯一比上一年下降的年份,所以也顯得異常,再有就是第33個數(shù)據(jù)(2012年)也顯示異常,說明去年的糧食產(chǎn)量數(shù)據(jù)較大,或者說距離均值較遠(yuǎn),也說明我國糧食產(chǎn)量持續(xù)增高。

3 結(jié)論

(1)在對糧食產(chǎn)量有眾多影響的因素中,經(jīng)過降維分析,最終我們選取出了兩個主要的因素,即x5：農(nóng)用化肥施用折純量(萬噸)和x7：糧食作物播種面積(千公頃),同時為了能夠比較這兩個因素對糧食產(chǎn)量的影響大小,我們對數(shù)據(jù)進(jìn)行了中心化處理。

(2)因為樣本觀測值只有33個,實際中又不可能增加,為了得到待估計系數(shù)的優(yōu)良統(tǒng)計性質(zhì),我們采用了近些年來較為流行的bootstrap方法來增加樣本量,同時采用較為穩(wěn)健的M-估計方法,得到了回歸結(jié)果。

(3)從回歸結(jié)果看,在固定一個自變量的情形下,農(nóng)用化肥施用折純量和糧食作物播種面積每增加一個單位,糧食產(chǎn)量將分別增加1.0645和0.3612個單位,說明化肥在糧食種植中的作用較大,在我國經(jīng)濟(jì)快速發(fā)展的形勢下,播種面積增加較為困難,重視化肥的使用將有助于糧食產(chǎn)量的增加。

(4)最后,基于Jackknife-after-Bootstrap方法,并結(jié)合殘差對糧食產(chǎn)量進(jìn)行了異常值分析,結(jié)果顯示觀測點(diǎn)1(1980年),2(1981年)增大了系數(shù) zx5,減小了系數(shù) zx7,而點(diǎn)33(2012年)則使兩個系數(shù)均增大。對照原始數(shù)據(jù)去看,因為我國糧食產(chǎn)量在近30年來增加較快,所以從數(shù)據(jù)上就使得1980年和1980年的糧食產(chǎn)量數(shù)據(jù)顯得異常,而第24個數(shù)據(jù)(2003年)是近些年來糧食產(chǎn)量唯一比上一年下降的年份,所以也顯得異常,再有就是第33個數(shù)據(jù)(2012年)也顯示異常,說明去年的糧食產(chǎn)量數(shù)據(jù)較大,或者說距離均值較遠(yuǎn),也說明我國糧食產(chǎn)量持續(xù)增高。

[1]Peter H,Joel H.A Simple Bootstrap Method for Constructing Nonparametric Confidence Bands for Functions[J].The Annals of Statistics,2013,41(4).

[2]Michael P.F,Erica H B,Michael A.P.Pointwise Confidence Intervals for A Survival Distribution With Small Samples or Heavy Censoring Biostat[J].Biostatistics,2013,14(4).

[3]Hoai T T,France M,Nicholas H.G.H.A Comparison of Bootstrap Approaches for Estimating Uncertainty of Parameters in Linear Mixed-Effects Models[J].Pharmaceutical Statistics,2013,12(3).

[4]劉薇，常振海.基于自助法的統(tǒng)計泛函估計比較研究[J].統(tǒng)計與決策,2013,(2).

[5]Huber P J.Robust Regression:Asymptotic Conjectures and Monte Carlo[J].Ann.Statist,1973,1(5).