r语言课程设计报告

-.z......资料....统计软件及应用课程设计虫情危害预测学院：班级：**：**：指导老师：目录一．背景与意义1二．问题重述1三．方法简介1(一).判别分析1(二).BP神经网络1四．数据处理与分析1(一).数据预处理1(二).判别分析11.模型建立12.模型优化13.模型应用1(三).BP神经网络11.模型建立12.模型优化13.建模重建14.模型应用1五．总结与建议1背景与意义农作物主要害虫常年对农作物造成严重危害，使农业经济遭到损失。预测害虫未来的发生动态，可以使治虫工作得以有目的、有计划、有重点的进行。害虫的预测预报工作是进行害虫综合防治的必要前提。只有对害虫发生危害的预测预报做到及时、准确，才能正确的拟定综合防治计划，及时采取必要的措施，经济有效的压低害虫的发生数量，保证农业的高产、稳产。问题重述本文选取的预测预报对象是庐江的田间水稻。水稻螟虫是水稻的重要害虫之一，对农作物的危害极大,其数量的多少一定程度上决定着水稻受危害的严重程度。通过对此昆虫的基本了解，发现气候因素对昆虫的发生发展有着密切关系，可以直接影响昆虫的生长、发育、生存、繁殖，从而造成害虫不同的发生期、发生量和危害程度。同时水稻螟虫是变温昆虫，其生长、发育和繁殖与气象条件的关系极为密切，所以我们从气候因素角度入手进行分析是合理的。本文从影响害虫生存繁殖的气候因素角度入手，结合往年的气象资料以及影响害虫生存繁殖的重要气候因素，选取了平均气温、最低气温、日照时间及降雨量四个主要影响因素，运用统计学方法确定虫害的发生量与气候因子的关系，并给出相应的预测方法。方法简介针对实际问题，我们需要通过对历史数据的分析，给出准则：当给定新时期下每一样本对应的各项气候指标时，能准确的判断其对应的虫害程度。下面介绍两种方法。(一).判别分析判别分析是在分类确定的条件下，根据*一研究对象的各种特征值判别其类型归属问题的一种多变量统计分析方法，是用以判别个体所属群体的一种统计方法。根据判别中的组数，可以分为两组判别分析和多组判别分析；根据判别函数的形式，可以分为线性判别和非线性判别；根据判别式处理变量的方法不同，可以分为逐步判别、序贯判别；根据判别标准不同，可以分为距离判别、Fisher判别、Bayes判别法。本文选用的是Fisher判别法。Fisher判别，是根据线性Fisher函数值进行判别，使用此准则要求各组变量的均值有显著性差异。该方法的基本思想是投影，即将原来在R维空间的自变量组合投影到维度较低的D维空间去，然后在D维空间中再进行分类。投影的原则是使得每一类的差异尽可能小，而不同类间投影的离差尽可能大。Fisher判别的优势在于对分布、方差等都没有任何限制，应用围比较广。另外，用该判别方法建立的判别方差可以直接用手工计算的方法进行新样品的判别，这在许多时候是非常方便的。(二).BP神经网络BP神经网络是1986年由Rumelhart和McCelland为首的科学家小组提出，是一种按误差逆传播算法训练的多层前馈网络，是目前应用最广泛的神经网络模型之一。BP网络能学习和存贮大量的输入-输出模式映射关系，而无需事前揭示描述这种映射关系的数学方程。它的学习规则是使用梯度下降法，通过反向传播来不断调整网络的权值和阈值，使网络的误差平方和最小。BP神经网络模型拓扑结构包括输入层、隐藏层和输出层。BP神经网络，由信息的正向传播和误差的反向传播两个过程组成。输入层各神经元负责接收来自外界的输入信息，并传递给中间层各神经元；中间层是部信息处理层，负责信息变换，根据信息变化能力的需求，中间层可以设计为单隐层或者多隐层结构；最后一个隐层传递到输出层各神经元的信息，经进一步处理后，完成一次学习的正向传播处理过程，由输出层向外界输出信息处理结果。当实际输出与期望输出不符时，进入误差的反向传播阶段。误差通过输出层，按误差梯度下降的方式修正各层权值，向隐层、输入层逐层反传。周而复始的信息正向传播和误差反向传播过程，是各层权值不断调整的过程，也是神经网络学习训练的过程，此过程一直进行到网络输出的误差减少到可以接受的程度，或者预先设定的学习次数为止。本文运用的是单层前馈网络模型，一般称为三层前馈网或三层感知器，即：输入层、中间层（也称隐层）和输出层。如下图所示它的特点是：各层神经元仅与相邻层神经元之间相互全连接，同层神经元之间无连接，各层神经元之间无反馈连接，构成具有层次结构的前馈型神经网络系统。数据处理与分析数据来自于1980～1988年6～10月庐江的月平均温度、月最低温度、月日照时间、月降雨量以及水稻的月平均虫害程度，且数据都是经过标准化之后的。(一).数据预处理首先，将降雨量的最后一个数据中的空格消除，再把原始数据的前两列合并为一列，并保存为CSV（逗号分隔）格式，之后用函数将数据读入到R的存中。>shuju<-read.csv('chongqing.csv')为了方便之后的操作，对各列进行命名，其中date代表年度时间，*1～*4分别代表平均气温、最低气温、日照时间、降雨量，y代表虫害程度。同时给出各类别下的样本量个数。由于本文不采用时间序列预测的方法进行分析，故第一列的数据使用不到，最终确定使用的数据集为剔除原数据框第一列之后的数据框。从原始数据中可知样本量为54，不是很多，如果再进行分块，可能由于样本量的限制无法得到最好的结果，因此，下文中的训练集和测试集均选用原始数据来充当。>names(shuju)<-c('date','*1','*2','*3','*4','y')>*<-shuju[,2:6]>table(y)y1234356211将数据框*中的变量到存中，降低后续代码的复杂程度，提高代码的可读性。>attach(*)(二).判别分析1.模型建立安装并加载软件包MASS。>install.packages(“MASS”)>library(MASS)将训练集中的y变量作为判别变量，其他剩余的变量作为特征变量，运用lda()函数来进行费希尔判别下的线性判别。>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输出结果的第二项为此次过程中各类别所使用的先验概率，其实就是训练集中判别变量各水平下的数量占总样本量的比例。第三项为各变量在每一类别中的均值，从中可以看出*1～*3在各类别下有较明显的差别，相反*4在各类别下的差别很小，也就说明降雨量对虫害程度的影响不是很显著。第四项给出线性判别式的参数矩阵。第五项给出各线性判别式分别的贡献比例。>plot(lda)从上图可以看到，在3个线性判别式下1～4这4个类别的分布情况，不同类别样本已经用相应数字标出。可能受样本量的限制，图中并不能看出很明显的分布情况。下面用上述模型对测试集进行预测，并给出测试集中y变量的预测结果与实际类别的混淆矩阵。>pred<-predict(lda,*[,1:4])>pred$class[1]11441111441111111111141111441111441111421[42]1114411111111Levels:1234>table(y,pred$class)y1234132003221033100146005上述混淆矩阵的行表示实际的类别，列表示预测判定的类别。在54个测试样本中，实际属于第1类的有35个，而由判定结果，在35个样本中，有32个判定正确，有3个被错判为第4类；第2类只有1个判定正确，有2个被错判为第1类，3个被错判为第4类；第3类没有一个判定正确；第4类有5个判定正确，6个被错判为第1类。从矩阵中只能看出每一类别的预测结果与实际类别的差距，只是一个直观的判断，并不能由此看出该模型的优劣，所以我们需要构造一个数量指标（误判概率=矩阵非对角线之和/样本总量），通过这个指标来确定模型的优劣。>e_lda<-sum(as.numeric(pred$class)!=y)/nrow(*)>e_lda[1]0.2962963通过计算，上述模型的误判率为29.6%，可以看出该模型的预测效果并不好，因此需要进一步的改进。2.模型优化从上文中我们知道，使用Fisher判别要求各组变量的均值有显著性差异，从上述结果已经知道降雨量在各类别下的均值差异性很小，所以我们可以考虑将其剔除，通过计算各变量在各类别下均值的方差，来作为剔除*一变量的理论依据。>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从结果来看，*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．模型应用如果运用上述两个模型对未来进行预测的话，只需将所需要的变量存放至一个新的数据框，需要注意的是，数据框中各变量的名字需要和建立模型时使用的数据框中各变量的名称一样，即平均气温、最低气温、日照时间、降雨量分别用*1～*4来表示。然后直接运用模型来给出预测结果。(三).BP神经网络1.模型建立安装并加载软件包nnet。>install.packages(“nnet”)>library(nnet)首先在建模之前，先介绍软件包nnet中的一个函数：class.ind()。该函数可以通过类别变量的因子变量来生成一个类指标矩阵。下面通过一个简单的例子来更清楚地表明该函数的功能。>z=c(1,2,3)*定义一个水平数为3的向量>class.ind(z)123[1,]100[2,]010[3,]001从输出结果可以看到，该函数主要是将向量变成一个矩阵，其中每行代表一个样本。只是将样本的类别用0和1来表示，即如果是该类，则在该类别名下用1表示，而其余的类别名下面用0表示。将数据框的前4列作为自变量数据，class.ind(y)作为类别变量数据，设置隐藏层节点个数为4，模型的最大迭代次数为1000次，用来防止模型的死循环。之后用核心函数nnet()来建立神经网络模型。>model=nnet(*[1:4],class.ind(y),data=*,size=4,+decay=5e-4,ma*it=1000)*神经网络>p=*[,1:4]*测试集定为原始数据>pred3=predict(model,p)*对测试集进行预测>head(pred3)*预测结果的部分显示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上述结果的每一行代表测试集中的每一个样本，其判别准则为：每一行中最大值所在的列名即为该样本所属的类别。>pnew=ma*.col(pred3)*确定测试集中每一个样本所属类别>table(y,pnew)*给出混淆矩阵pnewy12341331012240030011410010>e_e<-sum(pnew!=y)/nrow(*)*计算误判率>e_e[1]0.1111111从结果可以看出，该模型的误判率为11.1%，相比较判别分析大幅度降低，说明该模型对测试集的预测效果很好。>names(model)*模型中包含的所有输出结果[1]"n""nunits""nconn""conn""nsunits""decay"[7]"entropy""softma*""censored""value""wts""convergence"[13]"fitted.values""residuals""call"其中：value表示迭代最终值。wts表示模型在迭代过程中所寻找到的最优权重值，其取值为系统随机生成，即每次建模所使用的迭代初始值都是不相同的。由此可以看出神经网络模型的灵活性很高。convergence表示在模型建立的迭代过程中，迭代次数是否达到最大迭代次数。如果结果为1，则表明达到最大迭代次数；如果结果为0，则表明没有达到最大迭代次数。如果结果没有达到最大迭代次数，说明该模型可以被采纳，否则，就要增大参数ma*it的设定值，重新建模。>model$convergence[1]0该模型的convergence取值为0，说明该模型可以被采纳。2.模型优化虽然上述模型的各项指标可能都满足要求，但是，隐藏层中节点的个数和迭代次数都是人为设定的，随着两者数值的改变建立的模型也不尽相同，所以我们需要寻找两者的最优值，使其对应的模型最优。下面通过编写函数来实现这一目的。(1).寻找隐藏层节点个数最优值>f=function(n){+e=matri*(0,2,n,dimnames=list(c('误判概率',+'迭代终值'),c(1:n)))*定义一矩阵，用于存放结果+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='隐藏层结点个数',ylab+='误判率'))+}>f(10)*将n取为10输出结果：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结合矩阵以及折线图，可以看出当n=8的时候，误判率最小，同时迭代最终值也达到最小，说明当n=8即隐藏层节点个数为8时，模型的拟合效果最好，对测试集做出的预测效果也最好。(2).寻找迭代次数最优值>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='训练周+期',col='blue',ylim=c(min(e),ma*(e))))}ff(1000)从模型对测试集的误判率随训练周期的变化曲线可以看出，误判概率随着训练周期的增大先是急速下降，后趋于平稳。但是也有出现反向增长的现象，可能与训练集和测试集的选取有关。不过，从图中整体的变化趋势可以看出曲线的拐点在ma*it=200处，所以ma*it的最优值可能为200。3.建模重建取size=8，ma*it=200，重新建立神经网络模型，得到如下结果：>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%，较上次模型有所下降，但是模型的迭代次数达到了最大迭代次数，因此该模型可能还不是最优的。取size=8，ma*it=1000，重新建立神经网络模型，得到如下结果：>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(