r语言课程设计

1、统计软件及应用课程设计虫情危害预测学 院： 班 级： 学 号： 姓 名： 指导老师：目录一背景与意义 2.二问题重述 2.三方法简介 2.（ 一 ）. 判别分析 2.（ 二 ）.BP 神经网络 3.四数据处理与分析 3.（ 一 ）. 数据预处理 4.（ 二 ）. 判别分析 4.1. 模型建立 4.2. 模型优化 6.3. 模型应用 7.（ 三 ）.BP 神经网络 7.1. 模型建立 7.2. 模型优化 8.3. 建模重建 1.0.4. 模型应用 1.1.五总结与建议 1.2.一 背景与意义农作物主要害虫常年对农作物造成严重危害， 使农业经济遭到损失。 预测害 虫未来的发生动态，可以使治虫工作得

2、以有目的、有计划、有重点的进行。害虫 的预测预报工作是进行害虫综合防治的必要前提。 只有对害虫发生危害的预测预 报做到及时、准确，才能正确的拟定综合防治计划，及时采取必要的措施，经济 有效的压低害虫的发生数量，保证农业的高产、稳产。二 问题重述本文选取的预测预报对象是安徽庐江的田间水稻。 水稻螟虫是水稻的重要害 虫之一，对农作物的危害极大 ,其数量的多少一定程度上决定着水稻受危害的 严 重程度。通过对此昆虫的基本了解， 发现气候因素对昆虫的发生发展有着密切关 系，可以直接影响昆虫的生长、 发育、生存、繁殖，从而造成害虫不同的发生期、 发生量和危害程度。 同时水稻螟虫是变温昆虫， 其生长、 发育

3、和繁殖与气象条件 的关系极为密切，所以我们从气候因素角度入手进行分析是合理的。本文从影响害虫生存繁殖的气候因素角度入手， 结合往年的气象资料以及影 响害虫生存繁殖的重要气候因素， 选取了平均气温、 最低气温、 日照时间及降雨 量四个主要影响因素， 运用统计学方法确定虫害的发生量与气候因子的关系， 并 给出相应的预测方法。三 方法简介针对实际问题， 我们需要通过对历史数据的分析， 给出准则： 当给定新时期 下每一样本对应的各项气候指标时， 能准确的判断其对应的虫害程度。 下面介绍 两种方法。（ 一 ）. 判别分析判别分析是在分类确定的条件下， 根据某一研究对象的各种特征值判别其类 型归属问题的一

4、种多变量统计分析方法， 是用以判别个体所属群体的一种统计方 法。根据判别中的组数， 可以分为两组判别分析和多组判别分析； 根据判别函数 的形式，可以分为线性判别和非线性判别； 根据判别式处理变量的方法不同， 可 以分为逐步判别、序贯判别；根据判别标准不同，可以分为距离判别、 Fisher 判别、 Bayes 判别法。本文选用的是 Fisher 判别法。 Fisher 判别，是根据线性 Fisher 函数值进 行判别，使用此准则要求各组变量的均值有显著性差异。 该方法的基本思想是投 影，即将原来在 R维空间的自变量组合投影到维度较低的 D维空间去，然后在 D 维空间中再进行分类。 投影的原则是使

5、得每一类的差异尽可能小， 而不同类间投 影的离差尽可能大。 Fisher 判别的优势在于对分布、方差等都没有任何限制， 应用范围比较广。 另外，用该判别方法建立的判别方差可以直接用手工计算的方 法进行新样品的判别，这在许多时候是非常方便的。（ 二 ）.BP 神经网络BP神经网络是 1986年由 Rumelhart 和 McCelland 为首的科学家小组提出， 是一种按误差逆传播算法训练的多层前馈网络， 是目前应用最广泛的神经网络模 型之一。 BP 网络能学习和存贮大量的输入 - 输出模式映射关系，而无需事前揭示 描述这种映射关系的数学方程。 它的学习规则是使用梯度下降法， 通过反向传播 来不

6、断调整网络的权值和阈值，使网络的误差平方和最小。 BP神经网络模型拓 扑结构包括输入层、隐藏层和输出层。BP神经网络，由信息的正向传播和误差的反向传播两个过程组成。输入层 各神经元负责接收来自外界的输入信息， 并传递给中间层各神经元； 中间层是内 部信息处理层， 负责信息变换， 根据信息变化能力的需求， 中间层可以设计为单 隐层或者多隐层结构； 最后一个隐层传递到输出层各神经元的信息， 经进一步处 理后，完成一次学习的正向传播处理过程，由输出层向外界输出信息处理结果。 当实际输出与期望输出不符时， 进入误差的反向传播阶段。 误差通过输出层， 按 误差梯度下降的方式修正各层权值， 向隐层、 输入

7、层逐层反传。 周而复始的信息 正向传播和误差反向传播过程， 是各层权值不断调整的过程， 也是神经网络学习 训练的过程， 此过程一直进行到网络输出的误差减少到可以接受的程度， 或者预 先设定的学习次数为止。本文运用的是单层前馈网络模型，一般称为三层前馈网或三层感知器，即： 输入层、中间层（也称隐层）和输出层。如下图所示它的特点是：各层神经元仅 与相邻层神经元之间相互全连接， 同层内神经元之间无连接， 各层神经元之间无 反馈连接，构成具有层次结构的前馈型神经网络系统四 数据处理与分析数据来自于 19801988年 610 月安徽庐江的月平均温度、月最低温度、 月日照时间、 月降雨量以及水稻的月平均

8、虫害程度， 且数据都是经过标准化之后 的。( 一 ). 数据预处理首先，将降雨量的最后一个数据中的空格消除， 再把原始数据的前两列合并 为一列，并保存为 CSV(逗号分隔) 格式，之后用函数将数据读入到 R的内存中。 shuju names(shuju) x table(y)y1 2 3 435 6 2 11将数据框 X 中的变量链接到内存中， 降低后续代码的复杂程度， 提高代码的 可读性。 attach(x)( 二 ). 判别分析1. 模型建立安装并加载软件包 MASS。 install.packages(“ MASS”) library(MASS)将训练集中的 y 变量作为判别变量，其他剩

9、余的变量作为特征变量，运用 l da() 函数来进行费希尔判别下的线性判别。 lda ldaCall:lda(y ., data = x)Prior probabilities of groups:1 2 3 40.64814815 0.11111111 0.03703704 0.20370370Group 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

10、-0.49735004 0.56514545 0.608190909 0.2076455 -0.2936727Coefficients of linear discriminants: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 LD30.8767 0.1197 0.0036输出结果的第二项为此次过程

11、中各类别所使用的先验概率， 其实就是训练集 中判别变量各水平下的数量占总样本量的比例。 第三项为各变量在每一类别中的 均值，从中可以看出 x1x3在各类别下有较明显的差别， 相反 x4在各类别下的 差别很小， 也就说明降雨量对虫害程度的影响不是很显著。 第四项给出线性判别 式的参数矩阵。第五项给出各线性判别式分别的贡献比例。 plot(lda)从上图可以看到， 在 3 个线性判别式下 1 4 这 4 个类别的分布情况， 不同 类别样本已经用相应数字标出。 可能受样本量的限制， 图中并不能看出很明显的 分布情况。下面用上述模型对测试集进行预测， 并给出测试集中 y 变量的预测结果与实 际类别的混

12、淆矩阵。 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 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上述混淆矩阵的行表示实际的类别，列表示预测判定的类别。在 54 个测试 样本中，实际属于第 1 类的有 35 个，而由判定结果，在 35 个样本中，有 32 个判定正确，有 3

13、个被错判为第 4 类；第 2 类只有 1 个判定正确，有 2 个被错 判为第 1类，3个被错判为第 4类；第 3类没有一个判定正确；第 4 类有5个 判定正确， 6 个被错判为第 1 类。从矩阵中只能看出每一类别的预测结果与实际类别的差距， 只是一个直观的 判断，并不能由此看出该模型的优劣， 所以我们需要构造一个数量指标 (误判概 率 = 矩阵非对角线之和 / 样本总量)，通过这个指标来确定模型的优劣。 e_lda e_lda1 0.2962963通过计算，上述模型的误判率为 29.6% ，可以看出该模型的预测效果并不 好，因此需要进一步的改进。2. 模型优化从上文中我们知道，使用 Fishe

14、r 判别要求各组变量的均值有显著性差异， 从上述结果已经知道降雨量在各类别下的均值差异性很小， 所以我们可以考虑将 其剔除，通过计算各变量在各类别下均值的方差， 来作为剔除某一变量的理论依 据。 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从结果来看， x4 在各类别下均值的方差相对其他三者而言差别较大，因此 考虑将其剔除，然后重新进行线性判别分析。 lda2 pred2 table(y,pr

15、ed2$class)y 1 2 3 41 33 0 0 22 3 0 0 33 1 0 0 14 5 0 0 6 e_lda2 e_lda21 0.2777778从误判率的角度来看，剔除 x4 之后的模型误判率为 27.8%，相对原来的模 型有所降低；但从混淆矩阵来看，对于第 2类和第 3类的误判率为 100%，这相 对原来的模型更难让人接受。两个模型各有好坏，但整体而言，都不是很好。3模型应用如果运用上述两个模型对未来进行预测的话， 只需将所需要的变量存放至一 个新的数据框， 需要注意的是， 数据框中各变量的名字需要和建立模型时使用的 数据框中各变量的名称一样，即平均气温、最低气温、日照时间

16、、降雨量分别用 x1x4 来表示。然后直接运用模型来给出预测结果。( 三 ).BP 神经网络1. 模型建立安装并加载软件包 nnet 。 install.packages(“ nnet ”) library(nnet)首先在建模之前，先介绍软件包 nnet 中的一个函数： class.ind() 。该 函数可以通过类别变量的因子变量来生成一个类指标矩阵。 下面通过一个简单的 例子来更清楚地表明该函数的功能。 z=c(1,2,3) # 定义一个水平数为 3 的向量 class.ind(z)1 2 31, 1 0 02, 0 1 03, 0 0 1 从输出结果可以看到， 该函数主要是将向量变成一个

17、矩阵， 其中每行代表 个样本。只是将样本的类别用 0 和 1 来表示，即如果是该类，则在该类别名下 用 1 表示，而其余的类别名下面用 0 表示。将数据框的前 4 列作为自变量数据， class.ind(y) 作为类别变量数据， 设置隐藏层节点个数为 4 ，模型的最大迭代次数为 1000 次，用来防止模型的死 循环。之后用核心函数 nnet() 来建立神经网络模型。 model=nnet(x1:4,class.ind(y),data=x,size=4,# 神经网络 测试集定为原始数据 对测试集进行预测 预测结果的部分显示+ decay=5e-4,maxit=1000) p=x,1:4 # pr

18、ed3=predict(model,p) # head(pred3)#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, 0.8682759854 0.0000000000 5.810529e-03

19、1.975710e-016, 0.9999831991 0.0000000000 1.631513e-04 2.351911e-04上述结果的每一行代表测试集中的每一个样本， 其判别准则为： 每一行中最大值所在的列名即为该样本所属的类别。 pnew=max.col(pred3) #确定测试集中每一个样本所属类别 table(y,pnew)#给出混淆矩阵pnewy 1 2 3 41 33 1 0 12 2 4 0 03 0 0 1 14 1 0 0 10 e_e e_e1 0.1111111从结果可以看出，该模型的误判率为11.1% ，相比较判别分析大幅度降低，说明该模型对测试集的预测效果很好。

20、 names(model) # 模型中包含的所有输出结果1 n nunits nconn conn nsunits decay 7entropy softmax censored value wtsconvergence13 fitted.values residuals call其中： value 表示迭代最终值。wts 表示模型在迭代过程中所寻找到的最优权重值，其取值为系统随机生 成，即每次建模所使用的迭代初始值都是不相同的。 由此可以看出神经网络模型 的灵活性很高。convergence 表示在模型建立的迭代过程中， 迭代次数是否达到最大迭代 次数。如果结果为 1 ，则表明达到最大迭代次

21、数；如果结果为 0，则表明没有达 到最大迭代次数。 如果结果没有达到最大迭代次数， 说明该模型可以被采纳， 否 则，就要增大参数 maxit 的设定值，重新建模。 model$convergence1 0该模型的 convergence 取值为 0 ，说明该模型可以被采纳。2. 模型优化虽然上述模型的各项指标可能都满足要求， 但是，隐藏层中节点的个数和迭 代次数都是人为设定的， 随着两者数值的改变建立的模型也不尽相同， 所以我们 需要寻找两者的最优值， 使其对应的模型最优。 下面通过编写函数来实现这一目 的。(1). 寻找隐藏层节点个数最优值 f=function(n)+ e=matrix(0

22、,2,n,dimnames=list(c( 误判概率 ,+ 迭代终值 ),c(1:n) # 定义一矩阵，用于存放结果 + 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 f(10) # 将 n 取为 10 输出结果：1 2 3 4 5误判概率 0.2592593 0.2037037 0.1851852 0.0740741 0.0370371迭代终值 20

23、.707012 15.369598 15.779999 8.2334279 6.50628536 7 8 9 10误判概率 0.0185185 0.0370371 0.000000 0.0370371 0.0185185迭代终值 4.4975195 4.6649832 2.733296 4.3985989 3.2961255结合矩阵以及折线图， 可以看出当 n=8的时候，误判率最小， 同时迭代最终 值也达到最小， 说明当 n=8 即隐藏层节点个数为 8 时，模型的拟合效果最好， 对 测试集做出的预测效果也最好。(2). 寻找迭代次数最优值ff=function(n)+ e=0+ for(i i

24、n 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= 训练周+期 ,col=blue,ylim=c(min(e),max(e)ff(1000)从模型对测试集的误判率随训练周期的变化曲线可以看出， 误判概率随着训 练周期的增大先是急速下降， 后趋于平稳。 但是也有出现反向增长的现象， 可能 与训练集和测试集的选取有关。

25、 不过，从图中整体的变化趋势可以看出曲线的拐 点在 maxit=200 处，所以 maxit 的最优值可能为 200 。3. 建模重建取 size=8，maxit=200 ，重新建立神经网络模型，得到如下结果： 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.05555556 model1$convergence1 1误判概率为 5.6% ，较上次模型有所下降，但是模型的迭代次数达到了最大

26、迭代次数，因此该模型可能还不是最优的。取 size=8， maxit=1000，重新建立神经网络模型，得到如下结果： 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 e_e21 0 model2$convergence1 0从输出结果来看，该模型的误判率为 0，迭代次数也没有达到最大值，说明 该模型可以被采纳， 在样本量如此小的条件下， 建立的模型能有如此高的预测精 度，确实不得不让人佩服神经网络模型的强大。然而，由于训练集样本和测试集样本完全相同， 所以两者的特征也完全一样， 可能由于这个原因， 才会使得测试集的预测结果与实际值重合度如