r語言課程設(shè)計

1、統(tǒng)計軟件及應(yīng)用課程設(shè)計-一蟲情危害預(yù)測學 院：班 級：學 號：姓 名：指導老師:目錄1 .背景與意義2.2 .問題重述2.3 .方法簡介2.（1） .判別分析2.（2） .BP神經(jīng)網(wǎng)絡(luò)3.4 .數(shù)據(jù)處理與分析 3.（1） .數(shù)據(jù)預(yù)處理4.（2） .判別分析4.1 .模型建立4.2 .模型優(yōu)化6.3 .模型應(yīng)用7.（3） .BP神經(jīng)網(wǎng)絡(luò)7.1 .模型建立7.2 .模型優(yōu)化8.3 .建模重建1.04 .模型應(yīng)用1.15 .總結(jié)與建議1.2背景與意義農(nóng)作物主要害蟲常年對農(nóng)作物造成嚴重危害，使農(nóng)業(yè)經(jīng)濟遭到損失。預(yù)測害 蟲未來的發(fā)生動態(tài)，可以使治蟲工作得以有目的、有計劃、有重點的進行。害蟲 的預(yù)測預(yù)報工

2、作是進行害蟲綜合防治的必要前提。只有對害蟲發(fā)生危害的預(yù)測預(yù) 報做到及時、準確，才能正確的擬定綜合防治計劃，及時采取必要的措施，經(jīng)濟 有效的壓低害蟲的發(fā)生數(shù)量，保證農(nóng)業(yè)的高產(chǎn)、穩(wěn)產(chǎn)。2 .問題重述本文選取的預(yù)測預(yù)報對象是安徽廬江的田間水稻。 水稻螟蟲是水稻的重要害 蟲之一，對農(nóng)作物的危害極大，其數(shù)量的多少一定程度上決定著水稻受危害的嚴 重程度。通過對此昆蟲的基本了解，發(fā)現(xiàn)氣候因素對昆蟲的發(fā)生發(fā)展有著密切關(guān) 系，可以直接影響昆蟲的生長、發(fā)育、生存、繁殖，從而造成害蟲不同的發(fā)生期、 發(fā)生量和危害程度。同時水稻螟蟲是變溫昆蟲，其生長、發(fā)育和繁殖與氣象條件 的關(guān)系極為密切，所以我們從氣候因素角度入手進行

3、分析是合理的。本文從影響害蟲生存繁殖的氣候因素角度入手， 結(jié)合往年的氣象資料以及影 響害蟲生存繁殖的重要氣候因素， 選取了平均氣溫、最低氣溫、日照時間及降雨 量四個主要影響因素，運用統(tǒng)計學方法確定蟲害的發(fā)生量與氣候因子的關(guān)系，并給出相應(yīng)的預(yù)測方法。3 .方法簡介針對實際問題，我們需要通過對歷史數(shù)據(jù)的分析， 給出準則：當給定新時期 下每一樣本對應(yīng)的各項氣候指標時，能準確的判斷其對應(yīng)的蟲害程度。下面介紹 兩種方法。（1） .判別分析判別分析是在分類確定的條件下，根據(jù)某一研究對象的各種特征值判別其類 型歸屬問題的一種多變量統(tǒng)計分析方法，是用以判別個體所屬群體的一種統(tǒng)計方 法。根據(jù)判別中的組數(shù)，可以分

4、為兩組判別分析和多組判別分析；根據(jù)判別函數(shù)的形式，可以分為線性判別和非線性判別；根據(jù)判別式處理變量的方法不同， 可以分為逐步判別、序貫判別；根據(jù)判別標準不同，可以分為距離判別、Fisher判別、Bayes判別法。本文選用的是Fisher判別法。Fisher判別，是根據(jù)線性Fisher函數(shù)值進 行判別，使用此準則要求各組變量的均值有顯著性差異。該方法的基本思想是投 影，即將原來在R維空間的自變量組合投影到維度較低的 D維空間去，然后在D 維空間中再進行分類。投影的原則是使得每一類的差異盡可能小，而不同類間投 影的離差盡可能大。Fisher判別的優(yōu)勢在于對分布、方差等都沒有任何限制， 應(yīng)用范圍比較

5、廣。另外，用該判別方法建立的判別方差可以直接用手工計算的方 法進行新樣品的判別，這在許多時候是非常方便的。（2） .BP神經(jīng)網(wǎng)絡(luò)BP神經(jīng)網(wǎng)絡(luò)是1986年由Rumelhart和McCelland為首的科學家小組提出， 是一種按誤差逆?zhèn)鞑ニ惴ㄓ柧毜亩鄬忧梆伨W(wǎng)絡(luò)， 是目前應(yīng)用最廣泛的神經(jīng)網(wǎng)絡(luò)模 型之一。BP網(wǎng)絡(luò)能學習和存貯大量的輸入-輸出模式映射關(guān)系，而無需事前揭示 描述這種映射關(guān)系的數(shù)學方程。它的學習規(guī)則是使用梯度下降法，通過反向傳播 來不斷調(diào)整網(wǎng)絡(luò)的權(quán)值和閾值，使網(wǎng)絡(luò)的誤差平方和最小。BP神經(jīng)網(wǎng)絡(luò)模型拓撲結(jié)構(gòu)包括輸入層、隱藏層和輸出層。BP神經(jīng)網(wǎng)絡(luò)，由信息的正向傳播和誤差的反向傳播兩個過程組成。

6、輸入層 各神經(jīng)元負責接收來自外界的輸入信息， 并傳遞給中間層各神經(jīng)元；中間層是內(nèi) 部信息處理層，負責信息變換，根據(jù)信息變化能力的需求，中間層可以設(shè)計為單 隱層或者多隱層結(jié)構(gòu)；最后一個隱層傳遞到輸出層各神經(jīng)元的信息，經(jīng)進一步處理后，完成一次學習的正向傳播處理過程，由輸出層向外界輸出信息處理結(jié)果。 當實際輸出與期望輸出不符時，進入誤差的反向傳播階段。誤差通過輸出層，按 誤差梯度下降的方式修正各層權(quán)值，向隱層、輸入層逐層反傳。周而復(fù)始的信息 正向傳播和誤差反向傳播過程，是各層權(quán)值不斷調(diào)整的過程，也是神經(jīng)網(wǎng)絡(luò)學習 訓練的過程，此過程一直進行到網(wǎng)絡(luò)輸出的誤差減少到可以接受的程度，或者預(yù)先設(shè)定的學習次數(shù)為

7、止。本文運用的是單層前饋網(wǎng)絡(luò)模型，一般稱為三層前饋網(wǎng)或三層感知器，即： 輸入層、中間層（也稱隱層）和輸出層。如下圖所示它的特點是：各層神經(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ù)據(jù)都是經(jīng)過標準化之后 的。(1) .數(shù)據(jù)預(yù)處理首先，將降雨量的最后一個數(shù)據(jù)中的空格消除，再把原始數(shù)據(jù)的前兩列合并為一列，并保存為CSV逗號分隔)格式，之后用函數(shù)將數(shù)據(jù)讀入到 R的內(nèi)存中。>sh

8、uju<-read.csv('chongqing.csv')為了方便之后的操作，對各列進行命名，其中 date代表年度時間，x1x4 分別代表平均氣溫、最低氣溫、日照時間、降雨量，y代表蟲害程度。同時給出各類別下的樣本量個數(shù)。由于本文不采用時間序列預(yù)測的方法進行分析，故第一列的數(shù)據(jù)使用不到，最終確定使用的數(shù)據(jù)集為剔除原數(shù)據(jù)框第一列之后的數(shù)據(jù)框。 從原始數(shù)據(jù)中可知樣本量為54,不是很多，如果再進行分塊，可能由于樣本量 的限制無法得到最好的結(jié)果，因此，下文中的訓練集和測試集均選用原始數(shù)據(jù)來 充當。>names(shuju)<-c('date',&#

9、39;x1','x2','x3','x4','y')>x<-shuju,2:6>table(y)y1 2 3 435 6 2 11將數(shù)據(jù)框X中的變量鏈接到內(nèi)存中，降低后續(xù)代碼的復(fù)雜程度，提高代碼的 可讀性。>attach(x)(2) .判別分析1.模型建立安裝并加載軟件包MASS>install.packages( "MASS )>library(MASS)將訓練集中的y變量作為判別變量，其他剩余的變量作為特征變量，運用 l da()函數(shù)來進行費希爾判別下的線性判別。>l

10、da<-lda(y.,data=x) #線性判別>ldaCall:lda(y ., data = x)Prior probabilities of groups:12340.64814815 0.11111111 0.03703704 0.20370370Group means:x1 x2 x3 x41 -0.03336857-0.007231429-0.0008800-0.22869712 0.769233330.7488000000.1302167-0.26963333 0.377650000.4366000000.6406500-0.49735004 0.565145450.

11、6081909090.2076455-0.2936727Coefficients of linear discriminants:LD1 LD2 LD3 x1 1.1966105 -1.8868126 -3.1952390x2 0.9773588 1.1840427 3.3471367x3 0.2431136 2.1779657 -0.2219837x4 -0.2754396 -0.3782063 0.4283132Proportion of trace:LD1 LD2 LD30.8767 0.1197 0.0036輸出結(jié)果的第二項為此次過程中各類別所使用的先驗概率，其實就是訓練集中判別變量各

12、水平下的數(shù)量占總樣本量的比例。第三項為各變量在每一類別中的均值，從中可以看出x1x3在各類別下有較明顯的差別，相反x4在各類別下的 差別很小，也就說明降雨量對蟲害程度的影響不是很顯著。第四項給出線性判別 式的參數(shù)矩陣。第五項給出各線性判別式分別的貢獻比例。從上圖可以看到，在3個線性判別式下14這4個類別的分布情況，不同 類別樣本已經(jīng)用相應(yīng)數(shù)字標出?？赡苁軜颖玖康南拗?，圖中并不能看出很明顯的 分布情況。下面用上述模型對測試集進行預(yù)測，并給出測試集中y變量的預(yù)測結(jié)果與實 際類別的混淆矩陣。>pred<-predict(lda,x,1:4)>pred$class1 1 1 4 4

13、1 1 1 1 4 4 1 1 1 1 1 1 1 1 1 1 1 4 1 1 1 1 4 41 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 41 32 0 0 32 2 1 0 33 1 0 0 14 6 0 0 5上述混淆矩陣的行表示實際的類別，列表示預(yù)測判定的類別。在 54個測試 樣本中，實際屬于第1類的有35個，而由判定結(jié)果，在35個樣本中，有32 個判定正確，有3個被錯判為第4類；第2類只有1個判定正確，有2個被錯 判為第1類，3個被錯

14、判為第4類；第3類沒有一個判定正確；第4類有5個 判定正確，6個被錯判為第1類。從矩陣中只能看出每一類別的預(yù)測結(jié)果與實際類別的差距，只是一個直觀的判斷，并不能由此看出該模型的優(yōu)劣，所以我們需要構(gòu)造一個數(shù)量指標(誤判概 率二矩陣非對角線之和/樣本總量)，通過這個指標來確定模型的優(yōu)劣。>e_lda<-sum(as.numeric(pred$class)!=y)/nrow(x)>e_lda1 0.2962963通過計算，上述模型的誤判率為 29.6% ,可以看出該模型的預(yù)測效果并不 好，因此需要進一步的改進。2.模型優(yōu)化從上文中我們知道，使用Fisher判別要求各組變量的均值有顯著

15、性差異， 從上述結(jié)果已經(jīng)知道降雨量在各類別下的均值差異性很小，所以我們可以考慮將其剔除，通過計算各變量在各類別下均值的方差，來作為剔除某一變量的理論依據(jù)。>var=>c(var(lda$means,1),var(lda$means,2),+ var(lda$means,3),var(lda$means,4)>var1 0.11678909 0.10783314 0.07718836 0.01433243從結(jié)果來看，x4在各類別下均值的方差相對其他三者而言差別較大，因此 考慮將其剔除，然后重新進行線性判別分析。>lda2<-lda(yx1+x2+x3,data=x)

16、#線性判另 U>pred2<-predict(lda2,x,1:3)>table(y,pred2$class)y 1 2 3 41 33 0 0 22 3 0 0 33 1 0 0 14 5 0 0 6>e_lda2<-sum(as.numeric(pred2$class)!=y)/nrow(x)>e_lda21 0.2777778從誤判率的角度來看，剔除x4之后的模型誤判率為27.8%,相對原來的模 型有所降低；但從混淆矩陣來看，對于第 2類和第3類的誤判率為100%這相 對原來的模型更難讓人接受。兩個模型各有好壞，但整體而言，都不是很好。3.模型應(yīng)用如果

17、運用上述兩個模型對未來進行預(yù)測的話，只需將所需要的變量存放至一 個新的數(shù)據(jù)框，需要注意的是，數(shù)據(jù)框中各變量的名字需要和建立模型時使用的 數(shù)據(jù)框中各變量的名稱一樣，即平均氣溫、最低氣溫、日照時間、降雨量分別用 x1x4來表示。然后直接運用模型來給出預(yù)測結(jié)果。(三).BP神經(jīng)網(wǎng)絡(luò)1.模型建立安裝并加載軟件包nnet。>install.packages( “ nnet ” )>library(nnet)首先在建模之前，先介紹軟件包 nnet中的一個函數(shù)：class.ind() 。該 函數(shù)可以通過類別變量的因子變量來生成一個類指標矩陣。下面通過一個簡單的 例子來更清楚地表明該函數(shù)的功能。&

18、gt; z=c(1,2,3) #定義一個水平數(shù)為3的向量>class.ind(z)1 2 31, 1 0 02, 0 1 03, 0 0 1從輸出結(jié)果可以看到，該函數(shù)主要是將向量變成一個矩陣，其中每行代表一 個樣本。只是將樣本的類別用0和1來表示，即如果是該類，則在該類別名下用1表示，而其余的類別名下面用 0表示。將數(shù)據(jù)框的前4列作為自變量數(shù)據(jù)，class.ind(y)作為類別變量數(shù)據(jù)，設(shè)置隱藏層節(jié)點個數(shù)為4,模型的最大迭代次數(shù)為1000次，用來防止模型的死 循環(huán)。之后用核心函數(shù)nnet()來建立神經(jīng)網(wǎng)絡(luò)模型。>model=nnet(x1:4,class.ind(y),data=x

19、,size=4,+ decay=5e-4,maxit=1000) #神經(jīng)網(wǎng)絡(luò)> p=x,1:4#測試集定為原始數(shù)據(jù)> pred3=predict(model,p)#對測試集進行預(yù)測> head(pred3)#預(yù)測結(jié)果的部分顯示1 2341, 0.9999941595 0.0000000000 1.474405e-04 4.459695e-052, 0.9490308378 0.0113100101 1.454422e-05 9.221503e-033, 0.0049694310 0.0240625589 2.918010e-04 8.973849e-014, 0.000254

20、0397 0.0001941971 4.926029e-04 9.982468e-015, 0.8682759854 0.0000000000 5.810529e-03 1.975710e-016, 0.9999831991 0.0000000000 1.631513e-04 2.351911e-04上述結(jié)果的每一行代表測試集中的每一個樣本，其判別準則為：每一行中最大信所在的列名即為該樣本所屬的類別。>pnew=max.col(pred3) #確定測試集中每一個樣本所屬類別> table(y,pnew)#給出混淆矩陣pnewy 1 2 3 41 33 1 0 12 2 4 0 03

21、 0 0 1 14 1 0 0 10>e_e<-sum(pnew!=y)/nrow(x) #計算誤判率>e_e1 0.1111111從結(jié)果可以看出，該模型的誤判率為11.1% ,相比較判別分析大幅度降低， 說明該模型對測試集的預(yù)測效果很好。> names(model) #模型中包含的所有輸出結(jié)果1 "n" "nunits" "nconn" "conn" "nsunits" "decay" 7"entropy" "softm

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

23、次數(shù)；如果結(jié)果為0,則表明沒有達到最大迭代次數(shù)。如果結(jié)果沒有達到最大迭代次數(shù)，說明該模型可以被采納，否 則，就要增大參數(shù)maxit的設(shè)定值，重新建模。>model$convergence1 0該模型的convergence取值為0,說明該模型可以被采納。2 模型優(yōu)化雖然上述模型的各項指標可能都滿足要求，但是，隱藏層中節(jié)點的個數(shù)和迭代次數(shù)都是人為設(shè)定的，隨著兩者數(shù)值的改變建立的模型也不盡相同， 所以我們 需要尋找兩者的最優(yōu)值，使其對應(yīng)的模型最優(yōu)。下面通過編寫函數(shù)來實現(xiàn)這一目 的。(1) .尋找隱藏層節(jié)點個數(shù)最優(yōu)值> f=function(n)+ e=matrix(0,2,n,dimn

24、ames=list(c('誤判概率',+ 迭代終值'),c(1:n) # 定義一矩陣，用于存放結(jié)果+ for(i in 1: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<-sum(pnew!=y)/nrow(x)+e1,i=e_e+ e2,i=model$value+ + print(e)+ print(plot(as.numeric(e1,),

25、type='l',col='red',+ main='誤判率變化曲線,xlab='隱藏層結(jié)點個數(shù),ylab+=' 誤判率')+ > f(10)# 將 n 取為 10輸出結(jié)果：12345誤判概率0.2592593 0.2037037 0.18518520.0740741 0.0370371 迭代終值 20.707012 15.36959815.7799998.23342796.5062853678910誤判概率 0.01851850.03703710.000000 0.03703710.0185185迭代終值 4.497519

26、54.66498322.733296 4.39859893.2961255誤判率變化曲線結(jié)合矩陣以及折線圖，可以看出當n=8的時候，誤判率最小，同時迭代最終 值也達到最小，說明當n=8即隱藏層節(jié)點個數(shù)為8時，模型的擬合效果最好，對 測試集做出的預(yù)測效果也最好。(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

27、)+ +print(plot(1:n,e,T,ylab='誤判概率',xlab='訓練周+期',col='blue',ylim=c(min(e),max(e)ff(1000)woEd noO o020040060080010 00訓練周期從模型對測試集的誤判率隨訓練周期的變化曲線可以看出，誤判概率隨著訓練周期的增大先是急速下降，后趨于平穩(wěn)。但是也有出現(xiàn)反向增長的現(xiàn)象，可能 與訓練集和測試集的選取有關(guān)。不過，從圖中整體的變化趨勢可以看出曲線的拐 點在max讓=200 處，所以maxit的最優(yōu)值可能為200。(3) 模重建取size=8, maxit

28、=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<-sum(pnew1!=y)/nrow(x)> e_e11 0.05555556> model1$convergence1 1誤判概率為5.6% ,較上次模型有所下降，但是模型的迭代次數(shù)達到了最大 迭代次數(shù)，因此該模型可能還不是最優(yōu)的。取size=8, maxit

29、=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 41 35 0 0 02 0 6 0 03 0 0 2 04 0 0 0 11> e_e2<-sum(pnew2!=y)/nrow(x)> e_e21 0> model2$convergence1 0從輸出結(jié)果來看，該模型的誤判率為 0,迭代次數(shù)也沒有達到最大值，說明 該模型