r語言課程設(shè)計(jì)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上統(tǒng)計(jì)軟件及應(yīng)用課程設(shè)計(jì) -蟲情危害預(yù)測(cè) 學(xué) 院： 班 級(jí)： 學(xué) 號(hào)： 姓 名： 指導(dǎo)老師：目錄一 背景與意義農(nóng)作物主要害蟲常年對(duì)農(nóng)作物造成嚴(yán)重危害，使農(nóng)業(yè)經(jīng)濟(jì)遭到損失。預(yù)測(cè)害蟲未來的發(fā)生動(dòng)態(tài)，可以使治蟲工作得以有目的、有計(jì)劃、有重點(diǎn)的進(jìn)行。害蟲的預(yù)測(cè)預(yù)報(bào)工作是進(jìn)行害蟲綜合防治的必要前提。只有對(duì)害蟲發(fā)生危害的預(yù)測(cè)預(yù)報(bào)做到及時(shí)、準(zhǔn)確，才能正確的擬定綜合防治計(jì)劃，及時(shí)采取必要的措施，經(jīng)濟(jì)有效的壓低害蟲的發(fā)生數(shù)量，保證農(nóng)業(yè)的高產(chǎn)、穩(wěn)產(chǎn)。二 問題重述本文選取的預(yù)測(cè)預(yù)報(bào)對(duì)象是安徽廬江的田間水稻。水稻螟蟲是水稻的重要害蟲之一，對(duì)農(nóng)作物的危害極大,其數(shù)量的多少一定程度上決定著水稻受危

2、害的 嚴(yán)重程度。通過對(duì)此昆蟲的基本了解，發(fā)現(xiàn)氣候因素對(duì)昆蟲的發(fā)生發(fā)展有著密切關(guān)系，可以直接影響昆蟲的生長(zhǎng)、發(fā)育、生存、繁殖，從而造成害蟲不同的發(fā)生期、發(fā)生量和危害程度。同時(shí)水稻螟蟲是變溫昆蟲，其生長(zhǎng)、發(fā)育和繁殖與氣象條件的關(guān)系極為密切，所以我們從氣候因素角度入手進(jìn)行分析是合理的。本文從影響害蟲生存繁殖的氣候因素角度入手，結(jié)合往年的氣象資料以及影響害蟲生存繁殖的重要?dú)夂蛞蛩兀x取了平均氣溫、最低氣溫、日照時(shí)間及降雨量四個(gè)主要影響因素，運(yùn)用統(tǒng)計(jì)學(xué)方法確定蟲害的發(fā)生量與氣候因子的關(guān)系，并給出相應(yīng)的預(yù)測(cè)方法。三 方法簡(jiǎn)介針對(duì)實(shí)際問題，我們需要通過對(duì)歷史數(shù)據(jù)的分析，給出準(zhǔn)則：當(dāng)給定新時(shí)期下每一樣本對(duì)應(yīng)的

3、各項(xiàng)氣候指標(biāo)時(shí)，能準(zhǔn)確的判斷其對(duì)應(yīng)的蟲害程度。下面介紹兩種方法。(一).判別分析判別分析是在分類確定的條件下，根據(jù)某一研究對(duì)象的各種特征值判別其類型歸屬問題的一種多變量統(tǒng)計(jì)分析方法，是用以判別個(gè)體所屬群體的一種統(tǒng)計(jì)方法。根據(jù)判別中的組數(shù)，可以分為兩組判別分析和多組判別分析；根據(jù)判別函數(shù)的形式，可以分為線性判別和非線性判別；根據(jù)判別式處理變量的方法不同，可以分為逐步判別、序貫判別；根據(jù)判別標(biāo)準(zhǔn)不同，可以分為距離判別、Fisher判別、Bayes判別法。本文選用的是Fisher判別法。Fisher判別，是根據(jù)線性Fisher函數(shù)值進(jìn)行判別，使用此準(zhǔn)則要求各組變量的均值有顯著性差異。該方法的基本思想

4、是投影，即將原來在R維空間的自變量組合投影到維度較低的D維空間去，然后在D維空間中再進(jìn)行分類。投影的原則是使得每一類的差異盡可能小，而不同類間投影的離差盡可能大。Fisher判別的優(yōu)勢(shì)在于對(duì)分布、方差等都沒有任何限制，應(yīng)用范圍比較廣。另外，用該判別方法建立的判別方差可以直接用手工計(jì)算的方法進(jìn)行新樣品的判別，這在許多時(shí)候是非常方便的。(二).BP神經(jīng)網(wǎng)絡(luò)BP神經(jīng)網(wǎng)絡(luò)是1986年由Rumelhart和McCelland為首的科學(xué)家小組提出，是一種按誤差逆?zhèn)鞑ニ惴ㄓ?xùn)練的多層前饋網(wǎng)絡(luò)，是目前應(yīng)用最廣泛的神經(jīng)網(wǎng)絡(luò)模型之一。BP網(wǎng)絡(luò)能學(xué)習(xí)和存貯大量的輸入-輸出模式映射關(guān)系，而無需事前揭示描述這種映射關(guān)系的

5、數(shù)學(xué)方程。它的學(xué)習(xí)規(guī)則是使用梯度下降法，通過反向傳播來不斷調(diào)整網(wǎng)絡(luò)的權(quán)值和閾值，使網(wǎng)絡(luò)的誤差平方和最小。BP神經(jīng)網(wǎng)絡(luò)模型拓?fù)浣Y(jié)構(gòu)包括輸入層、隱藏層和輸出層。BP神經(jīng)網(wǎng)絡(luò)，由信息的正向傳播和誤差的反向傳播兩個(gè)過程組成。輸入層各神經(jīng)元負(fù)責(zé)接收來自外界的輸入信息，并傳遞給中間層各神經(jīng)元；中間層是內(nèi)部信息處理層，負(fù)責(zé)信息變換，根據(jù)信息變化能力的需求，中間層可以設(shè)計(jì)為單隱層或者多隱層結(jié)構(gòu)；最后一個(gè)隱層傳遞到輸出層各神經(jīng)元的信息，經(jīng)進(jìn)一步處理后，完成一次學(xué)習(xí)的正向傳播處理過程，由輸出層向外界輸出信息處理結(jié)果。當(dāng)實(shí)際輸出與期望輸出不符時(shí)，進(jìn)入誤差的反向傳播階段。誤差通過輸出層，按誤差梯度下降的方式修正各層權(quán)

6、值，向隱層、輸入層逐層反傳。周而復(fù)始的信息正向傳播和誤差反向傳播過程，是各層權(quán)值不斷調(diào)整的過程，也是神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)訓(xùn)練的過程，此過程一直進(jìn)行到網(wǎng)絡(luò)輸出的誤差減少到可以接受的程度，或者預(yù)先設(shè)定的學(xué)習(xí)次數(shù)為止。本文運(yùn)用的是單層前饋網(wǎng)絡(luò)模型，一般稱為三層前饋網(wǎng)或三層感知器，即：輸入層、中間層（也稱隱層）和輸出層。如下圖所示它的特點(diǎn)是：各層神經(jīng)元僅與相鄰層神經(jīng)元之間相互全連接，同層內(nèi)神經(jīng)元之間無連接，各層神經(jīng)元之間無反饋連接，構(gòu)成具有層次結(jié)構(gòu)的前饋型神經(jīng)網(wǎng)絡(luò)系統(tǒng)。四 數(shù)據(jù)處理與分析數(shù)據(jù)來自于19801988年610月安徽廬江的月平均溫度、月最低溫度、月日照時(shí)間、月降雨量以及水稻的月平均蟲害程度，且數(shù)據(jù)都

7、是經(jīng)過標(biāo)準(zhǔn)化之后的。(一).數(shù)據(jù)預(yù)處理首先，將降雨量的最后一個(gè)數(shù)據(jù)中的空格消除，再把原始數(shù)據(jù)的前兩列合并為一列，并保存為CSV（逗號(hào)分隔）格式，之后用函數(shù)將數(shù)據(jù)讀入到R的內(nèi)存中。 shuju names(shuju) x table(y)y 1 2 3 4 35 6 2 11 將數(shù)據(jù)框X中的變量鏈接到內(nèi)存中，降低后續(xù)代碼的復(fù)雜程度，提高代碼的可讀性。 attach(x)(二).判別分析1.模型建立安裝并加載軟件包MASS。 install.packages(“MASS”) library(MASS)將訓(xùn)練集中的y變量作為判別變量，其他剩余的變量作為特征變量，運(yùn)用lda()函數(shù)來進(jìn)行費(fèi)希爾判別下

8、的線性判別。 lda ldaCall:lda(y ., data = x)Prior probabilities of groups: 1 2 3 4 0. 0. 0. 0. Group means: x1 x2 x3 x41 -0. -0. -0. -0.2 0. 0. 0. -0.3 0. 0. 0. -0.4 0. 0. 0. -0.Coefficients of linear discriminants: LD1 LD2 LD3x1 1. -1. -3.x2 0. 1. 3.x3 0. 2. -0.x4 -0. -0. 0.Proportion of trace: LD1 LD2 LD

9、3 0.8767 0.1197 0.0036 輸出結(jié)果的第二項(xiàng)為此次過程中各類別所使用的先驗(yàn)概率，其實(shí)就是訓(xùn)練集中判別變量各水平下的數(shù)量占總樣本量的比例。第三項(xiàng)為各變量在每一類別中的均值，從中可以看出x1x3在各類別下有較明顯的差別，相反x4在各類別下的差別很小，也就說明降雨量對(duì)蟲害程度的影響不是很顯著。第四項(xiàng)給出線性判別式的參數(shù)矩陣。第五項(xiàng)給出各線性判別式分別的貢獻(xiàn)比例。 plot(lda)從上圖可以看到，在3個(gè)線性判別式下14這4個(gè)類別的分布情況，不同類別樣本已經(jīng)用相應(yīng)數(shù)字標(biāo)出?？赡苁軜颖玖康南拗?，圖中并不能看出很明顯的分布情況。下面用上述模型對(duì)測(cè)試集進(jìn)行預(yù)測(cè)，并給出測(cè)試集中y變量的預(yù)測(cè)結(jié)

10、果與實(shí)際類別的混淆矩陣。 pred pred$class1 1 1 4 4 1 1 1 1 4 4 1 1 1 1 1 1 1 1 1 1 1 4 1 1 1 1 4 4 1 1 1 1 4 4 1 1 1 1 4 2 142 1 1 1 4 4 1 1 1 1 1 1 1 1Levels: 1 2 3 4 table(y,pred$class)y 1 2 3 4 1 32 0 0 3 2 2 1 0 3 3 1 0 0 1 4 6 0 0 5上述混淆矩陣的行表示實(shí)際的類別，列表示預(yù)測(cè)判定的類別。在54個(gè)測(cè)試樣本中，實(shí)際屬于第1類的有35個(gè)，而由判定結(jié)果，在35個(gè)樣本中，有32個(gè)判定正確，有3

11、個(gè)被錯(cuò)判為第4類；第2類只有1個(gè)判定正確，有2個(gè)被錯(cuò)判為第1類，3個(gè)被錯(cuò)判為第4類；第3類沒有一個(gè)判定正確；第4類有5個(gè)判定正確，6個(gè)被錯(cuò)判為第1類。從矩陣中只能看出每一類別的預(yù)測(cè)結(jié)果與實(shí)際類別的差距，只是一個(gè)直觀的判斷，并不能由此看出該模型的優(yōu)劣，所以我們需要構(gòu)造一個(gè)數(shù)量指標(biāo)（誤判概率 = 矩陣非對(duì)角線之和/樣本總量），通過這個(gè)指標(biāo)來確定模型的優(yōu)劣。 e_lda e_lda1 0.通過計(jì)算，上述模型的誤判率為29.6%，可以看出該模型的預(yù)測(cè)效果并不好，因此需要進(jìn)一步的改進(jìn)。2.模型優(yōu)化從上文中我們知道，使用Fisher判別要求各組變量的均值有顯著性差異，從上述結(jié)果已經(jīng)知道降雨量在各類別下的均

12、值差異性很小，所以我們可以考慮將其剔除，通過計(jì)算各變量在各類別下均值的方差，來作為剔除某一變量的理論依據(jù)。 var=c(var(lda$means,1),var(lda$means,2),+ var(lda$means,3),var(lda$means,4) var1 0. 0. 0. 0.從結(jié)果來看，x4在各類別下均值的方差相對(duì)其他三者而言差別較大，因此考慮將其剔除，然后重新進(jìn)行線性判別分析。 lda2 pred2 table(y,pred2$class)y 1 2 3 4 1 33 0 0 2 2 3 0 0 3 3 1 0 0 1 4 5 0 0 6 e_lda2 e_lda21 0.從

13、誤判率的角度來看，剔除x4之后的模型誤判率為27.8%，相對(duì)原來的模型有所降低；但從混淆矩陣來看，對(duì)于第2類和第3類的誤判率為100%，這相對(duì)原來的模型更難讓人接受。兩個(gè)模型各有好壞，但整體而言，都不是很好。3模型應(yīng)用如果運(yùn)用上述兩個(gè)模型對(duì)未來進(jìn)行預(yù)測(cè)的話，只需將所需要的變量存放至一個(gè)新的數(shù)據(jù)框，需要注意的是，數(shù)據(jù)框中各變量的名字需要和建立模型時(shí)使用的數(shù)據(jù)框中各變量的名稱一樣，即平均氣溫、最低氣溫、日照時(shí)間、降雨量分別用x1x4來表示。然后直接運(yùn)用模型來給出預(yù)測(cè)結(jié)果。(三).BP神經(jīng)網(wǎng)絡(luò)1.模型建立安裝并加載軟件包nnet。 install.packages(“nnet”) library(n

14、net)首先在建模之前，先介紹軟件包nnet中的一個(gè)函數(shù)：class.ind()。該函數(shù)可以通過類別變量的因子變量來生成一個(gè)類指標(biāo)矩陣。下面通過一個(gè)簡(jiǎn)單的例子來更清楚地表明該函數(shù)的功能。 z=c(1,2,3) #定義一個(gè)水平數(shù)為3的向量 class.ind(z) 1 2 31, 1 0 02, 0 1 03, 0 0 1從輸出結(jié)果可以看到，該函數(shù)主要是將向量變成一個(gè)矩陣，其中每行代表一個(gè)樣本。只是將樣本的類別用0和1來表示，即如果是該類，則在該類別名下用1表示，而其余的類別名下面用0表示。將數(shù)據(jù)框的前4列作為自變量數(shù)據(jù)，class.ind(y)作為類別變量數(shù)據(jù)，設(shè)置隱藏層節(jié)點(diǎn)個(gè)數(shù)為4，模型的最

15、大迭代次數(shù)為1000次，用來防止模型的死循環(huán)。之后用核心函數(shù)nnet()來建立神經(jīng)網(wǎng)絡(luò)模型。 model=nnet(x1:4,class.ind(y),data=x,size=4,+ decay=5e-4,maxit=1000) #神經(jīng)網(wǎng)絡(luò) p=x,1:4 #測(cè)試集定為原始數(shù)據(jù) pred3=predict(model,p) #對(duì)測(cè)試集進(jìn)行預(yù)測(cè) head(pred3) #預(yù)測(cè)結(jié)果的部分顯示 1 2 3 41, 0. 0. 1.e-04 4.e-052, 0. 0. 1.e-05 9.e-033, 0. 0. 2.e-04 8.e-014, 0. 0. 4.e-04 9.e-015, 0. 0.

16、5.e-03 1.e-016, 0. 0. 1.e-04 2.e-04上述結(jié)果的每一行代表測(cè)試集中的每一個(gè)樣本，其判別準(zhǔn)則為：每一行中最大值所在的列名即為該樣本所屬的類別。 pnew=max.col(pred3) #確定測(cè)試集中每一個(gè)樣本所屬類別 table(y,pnew) #給出混淆矩陣 pnewy 1 2 3 4 1 33 1 0 1 2 2 4 0 0 3 0 0 1 1 4 1 0 0 10 e_e e_e1 0.從結(jié)果可以看出，該模型的誤判率為11.1%，相比較判別分析大幅度降低，說明該模型對(duì)測(cè)試集的預(yù)測(cè)效果很好。 names(model) #模型中包含的所有輸出結(jié)果1 n nuni

17、ts nconn conn nsunits decay 7 entropy softmax censored value wts convergence 13 fitted.values residuals call其中：value表示迭代最終值。wts表示模型在迭代過程中所尋找到的最優(yōu)權(quán)重值，其取值為系統(tǒng)隨機(jī)生成，即每次建模所使用的迭代初始值都是不相同的。由此可以看出神經(jīng)網(wǎng)絡(luò)模型的靈活性很高。convergence表示在模型建立的迭代過程中，迭代次數(shù)是否達(dá)到最大迭代次數(shù)。如果結(jié)果為1，則表明達(dá)到最大迭代次數(shù)；如果結(jié)果為0，則表明沒有達(dá)到最大迭代次數(shù)。如果結(jié)果沒有達(dá)到最大迭代次數(shù)，說明該模型可

18、以被采納，否則，就要增大參數(shù)maxit的設(shè)定值，重新建模。 model$convergence1 0該模型的convergence取值為0，說明該模型可以被采納。2.模型優(yōu)化雖然上述模型的各項(xiàng)指標(biāo)可能都滿足要求，但是，隱藏層中節(jié)點(diǎn)的個(gè)數(shù)和迭代次數(shù)都是人為設(shè)定的，隨著兩者數(shù)值的改變建立的模型也不盡相同，所以我們需要尋找兩者的最優(yōu)值，使其對(duì)應(yīng)的模型最優(yōu)。下面通過編寫函數(shù)來實(shí)現(xiàn)這一目的。(1).尋找隱藏層節(jié)點(diǎn)個(gè)數(shù)最優(yōu)值 f=function(n)+ e=matrix(0,2,n,dimnames=list(c(誤判概率,+ 迭代終值),c(1:n) #定義一矩陣，用于存放結(jié)果+ for(i in 1

19、:n)+ model=nnet(x1:4,class.ind(y),data=x,size=i,+ decay=5e-4,maxit=1000)+ pred=predict(model,p)+ pnew=max.col(pred)+ table(y,pnew)+ e_e f(10) #將n取為10輸出結(jié)果： 1 2 3 4 5 誤判概率 0. 0. 0. 0. 0.迭代終值 20. 15. 15. 8. 6. 6 7 8 9 10誤判概率 0. 0. 0. 0. 0.迭代終值 4. 4. 2. 4. 3.結(jié)合矩陣以及折線圖，可以看出當(dāng)n=8的時(shí)候，誤判率最小，同時(shí)迭代最終值也達(dá)到最小，說明當(dāng)n

20、=8即隱藏層節(jié)點(diǎn)個(gè)數(shù)為8時(shí)，模型的擬合效果最好，對(duì)測(cè)試集做出的預(yù)測(cè)效果也最好。(2).尋找迭代次數(shù)最優(yōu)值ff=function(n)+ e=0+ for(i in 1:n)+ model=nnet(x,1:4,class.ind(x$y),data=x,+ size=8,decay=5e-4,maxit=i)+ ei=sum(max.col(predict(model,x,1:4)!=+ x$y)/nrow(x)+ + print(plot(1:n,e,l,ylab=誤判概率,xlab=訓(xùn)練周+ 期,col=blue,ylim=c(min(e),max(e)ff(1000)從模型對(duì)測(cè)試集的誤判

21、率隨訓(xùn)練周期的變化曲線可以看出，誤判概率隨著訓(xùn)練周期的增大先是急速下降，后趨于平穩(wěn)。但是也有出現(xiàn)反向增長(zhǎng)的現(xiàn)象，可能與訓(xùn)練集和測(cè)試集的選取有關(guān)。不過，從圖中整體的變化趨勢(shì)可以看出曲線的拐點(diǎn)在maxit=200處，所以maxit的最優(yōu)值可能為200。3.建模重建取size=8，maxit=200，重新建立神經(jīng)網(wǎng)絡(luò)模型，得到如下結(jié)果： model1=nnet(x,1:4,class.ind(y),data=x,size=8,+ decay=5e-4,maxit=200) pred4=predict(model1,p) pnew1=max.col(pred4) e_e1 e_e11 0. model

22、1$convergence1 1誤判概率為5.6%，較上次模型有所下降，但是模型的迭代次數(shù)達(dá)到了最大迭代次數(shù)，因此該模型可能還不是最優(yōu)的。取size=8，maxit=1000，重新建立神經(jīng)網(wǎng)絡(luò)模型，得到如下結(jié)果： model2=nnet(x1:4,class.ind(y),data=x,size=8,+ decay=5e-4,maxit=1000) pred5=predict(model2,p) pnew2=max.col(pred5) table(y,pnew2) Pnew2y 1 2 3 4 1 35 0 0 0 2 0 6 0 0 3 0 0 2 0 4 0 0 0 11 e_e2 e_e21 0 model2$convergence1 0從輸出結(jié)果來看，該模型的誤判率為0，迭代次數(shù)也沒有達(dá)到最大值，說明該模型可以被采納，在樣本量如此小的條件下，建立的模型能有如此高的預(yù)測(cè)精度，確實(shí)不得不讓人佩服神經(jīng)網(wǎng)絡(luò)模型的強(qiáng)大。然而，由于訓(xùn)練集樣本和測(cè)試集樣本完全相同，所以兩者的特征也完全一樣，可能由于這個(gè)原因，才會(huì)使得測(cè)試集的預(yù)測(cè)結(jié)果與實(shí)際值重合度如此高；況且由于神經(jīng)網(wǎng)絡(luò)模型太過靈活，當(dāng)訓(xùn)練周期很長(zhǎng)的時(shí)候，建立出來的神經(jīng)網(wǎng)絡(luò)模型可能會(huì)記住訓(xùn)練集的幾乎所有信息，而測(cè)試集又和訓(xùn)練集相同，所以，模型的誤判率