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

1、統(tǒng)計(jì)軟件及應(yīng)用課程設(shè)計(jì) -蟲情危害預(yù)測 學(xué) 院： 班 級(jí)： 學(xué) 號(hào)： 姓 名： 指導(dǎo)老師：目錄一背景與意義2二問題重述2三方法簡介2(一).判別分析2(二).BP神經(jīng)網(wǎng)絡(luò)3四數(shù)據(jù)處理與分析3(一).數(shù)據(jù)預(yù)處理4(二).判別分析41.模型建立42.模型優(yōu)化63.模型應(yīng)用7(三).BP神經(jīng)網(wǎng)絡(luò)71.模型建立72.模型優(yōu)化83.建模重建104.模型應(yīng)用11五總結(jié)與建議12一 背景與意義農(nóng)作物主要害蟲常年對(duì)農(nóng)作物造成嚴(yán)重危害，使農(nóng)業(yè)經(jīng)濟(jì)遭到損失。預(yù)測害蟲未來的發(fā)生動(dòng)態(tài)，可以使治蟲工作得以有目的、有計(jì)劃、有重點(diǎn)的進(jìn)行。害蟲的預(yù)測預(yù)報(bào)工作是進(jìn)行害蟲綜合防治的必要前提。只有對(duì)害蟲發(fā)生危害的預(yù)測預(yù)報(bào)做到及時(shí)、

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

3、生存繁殖的重要?dú)夂蛞蛩?，選取了平均氣溫、最低氣溫、日照時(shí)間及降雨量四個(gè)主要影響因素，運(yùn)用統(tǒng)計(jì)學(xué)方法確定蟲害的發(fā)生量與氣候因子的關(guān)系，并給出相應(yīng)的預(yù)測方法。三 方法簡介針對(duì)實(shí)際問題，我們需要通過對(duì)歷史數(shù)據(jù)的分析，給出準(zhǔn)則：當(dāng)給定新時(shí)期下每一樣本對(duì)應(yīng)的各項(xiàng)氣候指標(biāo)時(shí)，能準(zhǔn)確的判斷其對(duì)應(yīng)的蟲害程度。下面介紹兩種方法。(一).判別分析判別分析是在分類確定的條件下，根據(jù)某一研究對(duì)象的各種特征值判別其類型歸屬問題的一種多變量統(tǒng)計(jì)分析方法，是用以判別個(gè)體所屬群體的一種統(tǒng)計(jì)方法。根據(jù)判別中的組數(shù)，可以分為兩組判別分析和多組判別分析；根據(jù)判別函數(shù)的形式，可以分為線性判別和非線性判別；根據(jù)判別式處理變量的方法不同

4、，可以分為逐步判別、序貫判別；根據(jù)判別標(biāo)準(zhǔn)不同，可以分為距離判別、Fisher判別、Bayes判別法。本文選用的是Fisher判別法。Fisher判別，是根據(jù)線性Fisher函數(shù)值進(jìn)行判別，使用此準(zhǔn)則要求各組變量的均值有顯著性差異。該方法的基本思想是投影，即將原來在R維空間的自變量組合投影到維度較低的D維空間去，然后在D維空間中再進(jìn)行分類。投影的原則是使得每一類的差異盡可能小，而不同類間投影的離差盡可能大。Fisher判別的優(yōu)勢(shì)在于對(duì)分布、方差等都沒有任何限制，應(yīng)用范圍比較廣。另外，用該判別方法建立的判別方差可以直接用手工計(jì)算的方法進(jìn)行新樣品的判別，這在許多時(shí)候是非常方便的。(二).BP神經(jīng)網(wǎng)

5、絡(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)系的數(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ù)信息變化能力的需求，中間層可

6、以設(shè)計(jì)為單隱層或者多隱層結(jié)構(gòu)；最后一個(gè)隱層傳遞到輸出層各神經(jīng)元的信息，經(jīng)進(jìn)一步處理后，完成一次學(xué)習(xí)的正向傳播處理過程，由輸出層向外界輸出信息處理結(jié)果。當(dāng)實(shí)際輸出與期望輸出不符時(shí)，進(jìn)入誤差的反向傳播階段。誤差通過輸出層，按誤差梯度下降的方式修正各層權(quán)值，向隱層、輸入層逐層反傳。周而復(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)

7、元之間相互全連接，同層內(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ù)都是經(jīng)過標(biāo)準(zhǔn)化之后的。(一).數(shù)據(jù)預(yù)處理首先，將降雨量的最后一個(gè)數(shù)據(jù)中的空格消除，再把原始數(shù)據(jù)的前兩列合并為一列，并保存為CSV（逗號(hào)分隔）格式，之后用函數(shù)將數(shù)據(jù)讀入到R的內(nèi)存中。> shuju<-read.csv('chongqing.csv') 為了方便之后的操作，對(duì)各列進(jìn)行命名，其中date代表年度時(shí)間，x1x4分別代

8、表平均氣溫、最低氣溫、日照時(shí)間、降雨量，y代表蟲害程度。同時(shí)給出各類別下的樣本量個(gè)數(shù)。由于本文不采用時(shí)間序列預(yù)測的方法進(jìn)行分析，故第一列的數(shù)據(jù)使用不到，最終確定使用的數(shù)據(jù)集為剔除原數(shù)據(jù)框第一列之后的數(shù)據(jù)框。從原始數(shù)據(jù)中可知樣本量為54，不是很多，如果再進(jìn)行分塊，可能由于樣本量的限制無法得到最好的結(jié)果，因此，下文中的訓(xùn)練集和測試集均選用原始數(shù)據(jù)來充當(dāng)。 > names(shuju)<-c('date','x1','x2','x3','x4','y')> x<-shuju,2:6

9、> 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)希爾判別下的線性判別。> lda<-lda(y.,data=x) #線性判別> ldaCall:lda(y ., data = x)Prior probabilities of g

10、roups: 1 2 3 4 0.64814815 0.11111111 0.03703704 0.20370370 Group means: x1 x2 x3 x41 -0.03336857 -0.007231429 -0.0008800 -0.22869712 0.76923333 0.748800000 0.1302167 -0.26963333 0.37765000 0.436600000 0.6406500 -0.49735004 0.56514545 0.608190909 0.2076455 -0.2936727Coefficients of linear discriminan

11、ts: LD1 LD2 LD3x1 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 LD3 0.8767 0.1197 0.0036 輸出結(jié)果的第二項(xiàng)為此次過程中各類別所使用的先驗(yàn)概率，其實(shí)就是訓(xùn)練集中判別變量各水平下的數(shù)量占總樣本量的比例。第三項(xiàng)為各變量在每一類別中的均值，從中可以看出x1x3在各類別下有較明顯的差別，相反x4

12、在各類別下的差別很小，也就說明降雨量對(duì)蟲害程度的影響不是很顯著。第四項(xiàng)給出線性判別式的參數(shù)矩陣。第五項(xiàng)給出各線性判別式分別的貢獻(xiàn)比例。> plot(lda)從上圖可以看到，在3個(gè)線性判別式下14這4個(gè)類別的分布情況，不同類別樣本已經(jīng)用相應(yīng)數(shù)字標(biāo)出?？赡苁軜颖玖康南拗疲瑘D中并不能看出很明顯的分布情況。下面用上述模型對(duì)測試集進(jìn)行預(yù)測，并給出測試集中y變量的預(yù)測結(jié)果與實(shí)際類別的混淆矩陣。> pred<-predict(lda,x,1:4)> 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

13、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ù)測判定的類別。在54個(gè)測試樣本中，實(shí)際屬于第1類的有35個(gè)，而由判定結(jié)果，在35個(gè)樣本中，有32個(gè)判定正確，有3個(gè)被錯(cuò)判為第4類；第2類只有1個(gè)判定正確，有2個(gè)被錯(cuò)判為第1類，3個(gè)被錯(cuò)判為第4類；第3類沒有一個(gè)判定正確；第4類有5個(gè)判定正確，6個(gè)被錯(cuò)判為第1類。從矩陣中只能看

14、出每一類別的預(yù)測結(jié)果與實(shí)際類別的差距，只是一個(gè)直觀的判斷，并不能由此看出該模型的優(yōu)劣，所以我們需要構(gòu)造一個(gè)數(shù)量指標(biāo)（誤判概率 = 矩陣非對(duì)角線之和/樣本總量），通過這個(gè)指標(biāo)來確定模型的優(yōu)劣。> e_lda<-sum(as.numeric(pred$class)!=y)/nrow(x)> e_lda1 0.2962963通過計(jì)算，上述模型的誤判率為29.6%，可以看出該模型的預(yù)測效果并不好，因此需要進(jìn)一步的改進(jìn)。2.模型優(yōu)化從上文中我們知道，使用Fisher判別要求各組變量的均值有顯著性差異，從上述結(jié)果已經(jīng)知道降雨量在各類別下的均值差異性很小，所以我們可以考慮將其剔除，通過計(jì)算

15、各變量在各類別下均值的方差，來作為剔除某一變量的理論依據(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在各類別下均值的方差相對(duì)其他三者而言差別較大，因此考慮將其剔除，然后重新進(jìn)行線性判別分析。> lda2<-lda(yx1+x2+x3,data=x) #線性判別> pred2<-predict(lda2,x,1:3)> t

16、able(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<-sum(as.numeric(pred2$class)!=y)/nrow(x)> e_lda21 0.2777778從誤判率的角度來看，剔除x4之后的模型誤判率為27.8%，相對(duì)原來的模型有所降低；但從混淆矩陣來看，對(duì)于第2類和第3類的誤判率為100%，這相對(duì)原來的模型更難讓人接受。兩個(gè)模型各有好壞，但整體而言，都不是很好。3模型應(yīng)用如果運(yùn)用上述兩個(gè)模型對(duì)未來進(jìn)行預(yù)測的話，只需將所需要的變量存放至一個(gè)新的數(shù)據(jù)框，需要注

17、意的是，數(shù)據(jù)框中各變量的名字需要和建立模型時(shí)使用的數(shù)據(jù)框中各變量的名稱一樣，即平均氣溫、最低氣溫、日照時(shí)間、降雨量分別用x1x4來表示。然后直接運(yùn)用模型來給出預(yù)測結(jié)果。(三).BP神經(jīng)網(wǎng)絡(luò)1.模型建立安裝并加載軟件包nnet。> install.packages(“nnet”)> library(nnet)首先在建模之前，先介紹軟件包nnet中的一個(gè)函數(shù)：class.ind()。該函數(shù)可以通過類別變量的因子變量來生成一個(gè)類指標(biāo)矩陣。下面通過一個(gè)簡單的例子來更清楚地表明該函數(shù)的功能。> z=c(1,2,3) #定義一個(gè)水平數(shù)為3的向量> class.ind(z) 1 2

18、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，模型的最大迭代次數(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

19、 #測試集定為原始數(shù)據(jù)> pred3=predict(model,p) #對(duì)測試集進(jìn)行預(yù)測> head(pred3) #預(yù)測結(jié)果的部分顯示 1 2 3 41, 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.0002540397 0.0001941971 4.926029e-04 9.982468e-015

20、, 0.8682759854 0.0000000000 5.810529e-03 1.975710e-016, 0.9999831991 0.0000000000 1.631513e-04 2.351911e-04上述結(jié)果的每一行代表測試集中的每一個(gè)樣本，其判別準(zhǔn)則為：每一行中最大值所在的列名即為該樣本所屬的類別。> pnew=max.col(pred3) #確定測試集中每一個(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<-sum(pn

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

22、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á)到最大

23、迭代次數(shù)，說明該模型可以被采納，否則，就要增大參數(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('誤判概率',+ '迭代終值

24、'),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,),type='l',col='red

25、9;,+ main='誤判率變化曲線',xlab='隱藏層結(jié)點(diǎn)個(gè)數(shù)',ylab+ ='誤判率')+ > f(10) #將n取為10輸出結(jié)果： 1 2 3 4 5 誤判概率 0.2592593 0.2037037 0.1851852 0.0740741 0.0370371迭代終值 20.707012 15.369598 15.779999 8.2334279 6.5062853 6 7 8 9 10誤判概率 0.0185185 0.0370371 0.000000 0.0370371 0.0185185迭代終值 4.4975195 4.664

26、9832 2.733296 4.3985989 3.2961255結(jié)合矩陣以及折線圖，可以看出當(dāng)n=8的時(shí)候，誤判率最小，同時(shí)迭代最終值也達(dá)到最小，說明當(dāng)n=8即隱藏層節(jié)點(diǎn)個(gè)數(shù)為8時(shí)，模型的擬合效果最好，對(duì)測試集做出的預(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)+ + print(

27、plot(1:n,e,'l',ylab='誤判概率',xlab='訓(xùn)練周+ 期',col='blue',ylim=c(min(e),max(e)ff(1000)從模型對(duì)測試集的誤判率隨訓(xùn)練周期的變化曲線可以看出，誤判概率隨著訓(xùn)練周期的增大先是急速下降，后趨于平穩(wěn)。但是也有出現(xiàn)反向增長的現(xiàn)象，可能與訓(xùn)練集和測試集的選取有關(guān)。不過，從圖中整體的變化趨勢(shì)可以看出曲線的拐點(diǎn)在maxit=200處，所以maxit的最優(yōu)值可能為200。3.建模重建取size=8，maxit=200，重新建立神經(jīng)網(wǎng)絡(luò)模型，得到如下結(jié)果：> model1=

28、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ù)達(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<-sum(pnew2!=y)/nrow(x)> e_e21 0> model2$convergence1 0從輸出結(jié)果來看，該模型的誤判率為0，迭代