r語言課程設(shè)計報告

-.z......資料....統(tǒng)計軟件及應(yīng)用課程設(shè)計蟲情危害預(yù)測學(xué)院：班級：**：**：指導(dǎo)老師：目錄一．背景與意義1二．問題重述1三．方法簡介1(一).判別分析1(二).BP神經(jīng)網(wǎng)絡(luò)1四．數(shù)據(jù)處理與分析1(一).數(shù)據(jù)預(yù)處理1(二).判別分析11.模型建立12.模型優(yōu)化13.模型應(yīng)用1(三).BP神經(jīng)網(wǎng)絡(luò)11.模型建立12.模型優(yōu)化13.建模重建14.模型應(yīng)用1五．總結(jié)與建議1背景與意義農(nóng)作物主要害蟲常年對農(nóng)作物造成嚴重危害，使農(nóng)業(yè)經(jīng)濟遭到損失。預(yù)測害蟲未來的發(fā)生動態(tài)，可以使治蟲工作得以有目的、有計劃、有重點的進行。害蟲的預(yù)測預(yù)報工作是進行害蟲綜合防治的必要前提。只有對害蟲發(fā)生危害的預(yù)測預(yù)報做到及時、準確，才能正確的擬定綜合防治計劃，及時采取必要的措施，經(jīng)濟有效的壓低害蟲的發(fā)生數(shù)量，保證農(nóng)業(yè)的高產(chǎn)、穩(wěn)產(chǎn)。問題重述本文選取的預(yù)測預(yù)報對象是廬江的田間水稻。水稻螟蟲是水稻的重要害蟲之一，對農(nóng)作物的危害極大,其數(shù)量的多少一定程度上決定著水稻受危害的嚴重程度。通過對此昆蟲的基本了解，發(fā)現(xiàn)氣候因素對昆蟲的發(fā)生發(fā)展有著密切關(guān)系，可以直接影響昆蟲的生長、發(fā)育、生存、繁殖，從而造成害蟲不同的發(fā)生期、發(fā)生量和危害程度。同時水稻螟蟲是變溫昆蟲，其生長、發(fā)育和繁殖與氣象條件的關(guān)系極為密切，所以我們從氣候因素角度入手進行分析是合理的。本文從影響害蟲生存繁殖的氣候因素角度入手，結(jié)合往年的氣象資料以及影響害蟲生存繁殖的重要氣候因素，選取了平均氣溫、最低氣溫、日照時間及降雨量四個主要影響因素，運用統(tǒng)計學(xué)方法確定蟲害的發(fā)生量與氣候因子的關(guān)系，并給出相應(yīng)的預(yù)測方法。方法簡介針對實際問題，我們需要通過對歷史數(shù)據(jù)的分析，給出準則：當給定新時期下每一樣本對應(yīng)的各項氣候指標時，能準確的判斷其對應(yīng)的蟲害程度。下面介紹兩種方法。(一).判別分析判別分析是在分類確定的條件下，根據(jù)*一研究對象的各種特征值判別其類型歸屬問題的一種多變量統(tǒng)計分析方法，是用以判別個體所屬群體的一種統(tǒng)計方法。根據(jù)判別中的組數(shù)，可以分為兩組判別分析和多組判別分析；根據(jù)判別函數(shù)的形式，可以分為線性判別和非線性判別；根據(jù)判別式處理變量的方法不同，可以分為逐步判別、序貫判別；根據(jù)判別標準不同，可以分為距離判別、Fisher判別、Bayes判別法。本文選用的是Fisher判別法。Fisher判別，是根據(jù)線性Fisher函數(shù)值進行判別，使用此準則要求各組變量的均值有顯著性差異。該方法的基本思想是投影，即將原來在R維空間的自變量組合投影到維度較低的D維空間去，然后在D維空間中再進行分類。投影的原則是使得每一類的差異盡可能小，而不同類間投影的離差盡可能大。Fisher判別的優(yōu)勢在于對分布、方差等都沒有任何限制，應(yīng)用圍比較廣。另外，用該判別方法建立的判別方差可以直接用手工計算的方法進行新樣品的判別，這在許多時候是非常方便的。(二).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)系的數(shù)學(xué)方程。它的學(xué)習(xí)規(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ò)，由信息的正向傳播和誤差的反向傳播兩個過程組成。輸入層各神經(jīng)元負責(zé)接收來自外界的輸入信息，并傳遞給中間層各神經(jīng)元；中間層是部信息處理層，負責(zé)信息變換，根據(jù)信息變化能力的需求，中間層可以設(shè)計為單隱層或者多隱層結(jié)構(gòu)；最后一個隱層傳遞到輸出層各神經(jīng)元的信息，經(jīng)進一步處理后，完成一次學(xué)習(xí)的正向傳播處理過程，由輸出層向外界輸出信息處理結(jié)果。當實際輸出與期望輸出不符時，進入誤差的反向傳播階段。誤差通過輸出層，按誤差梯度下降的方式修正各層權(quán)值，向隱層、輸入層逐層反傳。周而復(fù)始的信息正向傳播和誤差反向傳播過程，是各層權(quán)值不斷調(diào)整的過程，也是神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)訓(xùn)練的過程，此過程一直進行到網(wǎng)絡(luò)輸出的誤差減少到可以接受的程度，或者預(yù)先設(shè)定的學(xué)習(xí)次數(shù)為止。本文運用的是單層前饋網(wǎng)絡(luò)模型，一般稱為三層前饋網(wǎng)或三層感知器，即：輸入層、中間層（也稱隱層）和輸出層。如下圖所示它的特點是：各層神經(jīng)元僅與相鄰層神經(jīng)元之間相互全連接，同層神經(jīng)元之間無連接，各層神經(jīng)元之間無反饋連接，構(gòu)成具有層次結(jié)構(gòu)的前饋型神經(jīng)網(wǎng)絡(luò)系統(tǒng)。數(shù)據(jù)處理與分析數(shù)據(jù)來自于1980～1988年6～10月廬江的月平均溫度、月最低溫度、月日照時間、月降雨量以及水稻的月平均蟲害程度，且數(shù)據(jù)都是經(jīng)過標準化之后的。(一).數(shù)據(jù)預(yù)處理首先，將降雨量的最后一個數(shù)據(jù)中的空格消除，再把原始數(shù)據(jù)的前兩列合并為一列，并保存為CSV（逗號分隔）格式，之后用函數(shù)將數(shù)據(jù)讀入到R的存中。>shuju<-read.csv('chongqing.csv')為了方便之后的操作，對各列進行命名，其中date代表年度時間，*1～*4分別代表平均氣溫、最低氣溫、日照時間、降雨量，y代表蟲害程度。同時給出各類別下的樣本量個數(shù)。由于本文不采用時間序列預(yù)測的方法進行分析，故第一列的數(shù)據(jù)使用不到，最終確定使用的數(shù)據(jù)集為剔除原數(shù)據(jù)框第一列之后的數(shù)據(jù)框。從原始數(shù)據(jù)中可知樣本量為54，不是很多，如果再進行分塊，可能由于樣本量的限制無法得到最好的結(jié)果，因此，下文中的訓(xùn)練集和測試集均選用原始數(shù)據(jù)來充當。>names(shuju)<-c('date','*1','*2','*3','*4','y')>*<-shuju[,2:6]>table(y)y1234356211將數(shù)據(jù)框*中的變量到存中，降低后續(xù)代碼的復(fù)雜程度，提高代碼的可讀性。>attach(*)(二).判別分析1.模型建立安裝并加載軟件包MASS。>install.packages(“MASS”)>library(MASS)將訓(xùn)練集中的y變量作為判別變量，其他剩余的變量作為特征變量，運用lda()函數(shù)來進行費希爾判別下的線性判別。>lda<-lda(y~.,data=*)*線性判別>ldaCall:lda(y~.,data=*)Priorprobabilitiesofgroups:12340.648148150.111111110.037037040.20370370Groupmeans:*1*2*3*41-0.03336857-0.007231429-0.0008800-0.228697120.769233330.7488000000.1302167-0.269633330.377650000.4366000000.6406500-0.497350040.565145450.6081909090.2076455-0.2936727Coefficientsoflineardiscriminants:LD1LD2LD3*11.1966105-1.8868126-3.1952390*20.97735881.18404273.3471367*30.24311362.1779657-0.2219837*4-0.2754396-0.37820630.4283132Proportionoftrace:LD1LD2LD30.87670.11970.0036輸出結(jié)果的第二項為此次過程中各類別所使用的先驗概率，其實就是訓(xùn)練集中判別變量各水平下的數(shù)量占總樣本量的比例。第三項為各變量在每一類別中的均值，從中可以看出*1～*3在各類別下有較明顯的差別，相反*4在各類別下的差別很小，也就說明降雨量對蟲害程度的影響不是很顯著。第四項給出線性判別式的參數(shù)矩陣。第五項給出各線性判別式分別的貢獻比例。>plot(lda)從上圖可以看到，在3個線性判別式下1～4這4個類別的分布情況，不同類別樣本已經(jīng)用相應(yīng)數(shù)字標出?？赡苁軜颖玖康南拗疲瑘D中并不能看出很明顯的分布情況。下面用上述模型對測試集進行預(yù)測，并給出測試集中y變量的預(yù)測結(jié)果與實際類別的混淆矩陣。>pred<-predict(lda,*[,1:4])>pred$class[1]11441111441111111111141111441111441111421[42]1114411111111Levels:1234>table(y,pred$class)y1234132003221033100146005上述混淆矩陣的行表示實際的類別，列表示預(yù)測判定的類別。在54個測試樣本中，實際屬于第1類的有35個，而由判定結(jié)果，在35個樣本中，有32個判定正確，有3個被錯判為第4類；第2類只有1個判定正確，有2個被錯判為第1類，3個被錯判為第4類；第3類沒有一個判定正確；第4類有5個判定正確，6個被錯判為第1類。從矩陣中只能看出每一類別的預(yù)測結(jié)果與實際類別的差距，只是一個直觀的判斷，并不能由此看出該模型的優(yōu)劣，所以我們需要構(gòu)造一個數(shù)量指標（誤判概率=矩陣非對角線之和/樣本總量），通過這個指標來確定模型的優(yōu)劣。>e_lda<-sum(as.numeric(pred$class)!=y)/nrow(*)>e_lda[1]0.2962963通過計算，上述模型的誤判率為29.6%，可以看出該模型的預(yù)測效果并不好，因此需要進一步的改進。2.模型優(yōu)化從上文中我們知道，使用Fisher判別要求各組變量的均值有顯著性差異，從上述結(jié)果已經(jīng)知道降雨量在各類別下的均值差異性很小，所以我們可以考慮將其剔除，通過計算各變量在各類別下均值的方差，來作為剔除*一變量的理論依據(jù)。>var=>c(var(lda$means[,1]),var(lda$means[,2]),+var(lda$means[,3]),var(lda$means[,4]))>var[1]0.116789090.107833140.077188360.01433243從結(jié)果來看，*4在各類別下均值的方差相對其他三者而言差別較大，因此考慮將其剔除，然后重新進行線性判別分析。>lda2<-lda(y~*1+*2+*3,data=*)*線性判別>pred2<-predict(lda2,*[,1:3])>table(y,pred2$class)y1234133002230033100145006>e_lda2<-sum(as.numeric(pred2$class)!=y)/nrow(*)>e_lda2[1]0.2777778從誤判率的角度來看，剔除*4之后的模型誤判率為27.8%，相對原來的模型有所降低；但從混淆矩陣來看，對于第2類和第3類的誤判率為100%，這相對原來的模型更難讓人接受。兩個模型各有好壞，但整體而言，都不是很好。3．模型應(yīng)用如果運用上述兩個模型對未來進行預(yù)測的話，只需將所需要的變量存放至一個新的數(shù)據(jù)框，需要注意的是，數(shù)據(jù)框中各變量的名字需要和建立模型時使用的數(shù)據(jù)框中各變量的名稱一樣，即平均氣溫、最低氣溫、日照時間、降雨量分別用*1～*4來表示。然后直接運用模型來給出預(yù)測結(jié)果。(三).BP神經(jīng)網(wǎng)絡(luò)1.模型建立安裝并加載軟件包nnet。>install.packages(“nnet”)>library(nnet)首先在建模之前，先介紹軟件包nnet中的一個函數(shù)：class.ind()。該函數(shù)可以通過類別變量的因子變量來生成一個類指標矩陣。下面通過一個簡單的例子來更清楚地表明該函數(shù)的功能。>z=c(1,2,3)*定義一個水平數(shù)為3的向量>class.ind(z)123[1,]100[2,]010[3,]001從輸出結(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(*[1:4],class.ind(y),data=*,size=4,+decay=5e-4,ma*it=1000)*神經(jīng)網(wǎng)絡(luò)>p=*[,1:4]*測試集定為原始數(shù)據(jù)>pred3=predict(model,p)*對測試集進行預(yù)測>head(pred3)*預(yù)測結(jié)果的部分顯示1234[1,]0.99999415950.00000000001.474405e-044.459695e-05[2,]0.94903083780.01131001011.454422e-059.221503e-03[3,]0.00496943100.02406255892.918010e-048.973849e-01[4,]0.00025403970.00019419714.926029e-049.982468e-01[5,]0.86827598540.00000000005.810529e-031.975710e-01[6,]0.99998319910.00000000001.631513e-042.351911e-04上述結(jié)果的每一行代表測試集中的每一個樣本，其判別準則為：每一行中最大值所在的列名即為該樣本所屬的類別。>pnew=ma*.col(pred3)*確定測試集中每一個樣本所屬類別>table(y,pnew)*給出混淆矩陣pnewy12341331012240030011410010>e_e<-sum(pnew!=y)/nrow(*)*計算誤判率>e_e[1]0.1111111從結(jié)果可以看出，該模型的誤判率為11.1%，相比較判別分析大幅度降低，說明該模型對測試集的預(yù)測效果很好。>names(model)*模型中包含的所有輸出結(jié)果[1]"n""nunits""nconn""conn""nsunits""decay"[7]"entropy""softma*""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，則表明達到最大迭代次數(shù)；如果結(jié)果為0，則表明沒有達到最大迭代次數(shù)。如果結(jié)果沒有達到最大迭代次數(shù)，說明該模型可以被采納，否則，就要增大參數(shù)ma*it的設(shè)定值，重新建模。>model$convergence[1]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=matri*(0,2,n,dimnames=list(c('誤判概率',+'迭代終值'),c(1:n)))*定義一矩陣，用于存放結(jié)果+for(iin1:n){+model=nnet(*[1:4],class.ind(y),data=*,size=i,+decay=5e-4,ma*it=1000)+pred=predict(model,p)+pnew=ma*.col(pred)+table(y,pnew)+e_e<-sum(pnew!=y)/nrow(*)+e[1,i]=e_e+e[2,i]=model$value+}+print(e)+print(plot(as.numeric(e[1,]),type='l',col='red',+main='誤判率變化曲線',*lab='隱藏層結(jié)點個數(shù)',ylab+='誤判率'))+}>f(10)*將n取為10輸出結(jié)果：12345誤判概率0.25925930.20370370.18518520.07407410.0370371迭代終值20.70701215.36959815.7799998.23342796.5062853678910誤判概率0.01851850.03703710.0000000.03703710.0185185迭代終值4.49751954.66498322.7332964.39859893.2961255結(jié)合矩陣以及折線圖，可以看出當n=8的時候，誤判率最小，同時迭代最終值也達到最小，說明當n=8即隱藏層節(jié)點個數(shù)為8時，模型的擬合效果最好，對測試集做出的預(yù)測效果也最好。(2).尋找迭代次數(shù)最優(yōu)值>ff=function(n){+e=0+for(iin1:n){+model=nnet(*[,1:4],class.ind(*$y),data=*,+size=8,decay=5e-4,ma*it=i)+e[i]=sum(ma*.col(predict(model,*[,1:4]))!=+*$y)/nrow(*)+}+print(plot(1:n,e,'l',ylab='誤判概率',*lab='訓(xùn)練周+期',col='blue',ylim=c(min(e),ma*(e))))}ff(1000)從模型對測試集的誤判率隨訓(xùn)練周期的變化曲線可以看出，誤判概率隨著訓(xùn)練周期的增大先是急速下降，后趨于平穩(wěn)。但是也有出現(xiàn)反向增長的現(xiàn)象，可能與訓(xùn)練集和測試集的選取有關(guān)。不過，從圖中整體的變化趨勢可以看出曲線的拐點在ma*it=200處，所以ma*it的最優(yōu)值可能為200。3.建模重建取size=8，ma*it=200，重新建立神經(jīng)網(wǎng)絡(luò)模型，得到如下結(jié)果：>model1=nnet(*[,1:4],class.ind(y),data=*,size=8,+decay=5e-4,ma*it=200)>pred4=predict(model1,p)>pnew1=ma*.col(pred4)>e_e1<-sum(pnew1!=y)/nrow(*)>e_e1[1]0.05555556>model1$convergence[1]1誤判概率為5.6%，較上次模型有所下降，但是模型的迭代次數(shù)達到了最大迭代次數(shù)，因此該模型可能還不是最優(yōu)的。取size=8，ma*it=1000，重新建立神經(jīng)網(wǎng)絡(luò)模型，得到如下結(jié)果：>model2=nnet(*[1:4],class.ind(y),data=*,size=8,+decay=5e-4,ma*it=1000)>pred5=predict(model2,p)>pnew2=ma*.col(pred5)>table(y,pnew2)Pnew2y12341350002060030020400011>e_e2<-sum(pnew2!=y)/nrow(