统计软件与应用课程设计

统计软件与应用课程设计

----虫情危害预料

学院:

班级:

学号:

姓名:

指导老师:

书目

一.背景与意义....................................................2

问题重述......................................................2

三.方法简介.....................................................2

（一）.判别分析.................................................3

（二）.BP神经网络..............................................4

四.数据处理与分析...............................................5

（一）.数据预处理...............................................5

（二）.判别分析.................................................6

1.模型建立................................................6

2.模型优化...............................................10

3.模型应月...............................................11

（三）.BP神经网络.............................................11

1.模型建立...............................................11

2.模型优化...............................................14

3.建模重建...............................................18

4.模型应月...............................................20

五.总结与建议..................................................20

背景与意义

农作物主要害虫常年对农作物造成严峻危害，使农业经济遭到损失。

预料害虫将来的发生动态，可以使治虫工作得以有目的、有安排、有重点

的进行。害虫的预料预报工作是进行害虫综合防治的必要前提。只有对害

虫发生危害的预料预报做到与时、精确，才能正确的拟定综合防治安排,

与时实行必要的措施，经济有效的压低害虫的发生数量，保证农业的高产、

稳产。

问题重述

本文选取的预料预报对象是安徽庐江的田间水稻。水稻螟虫是水稻的

重要害虫之一，对农作物的危害极大，其数量的多少确定程度上确定着水稻

受危害的严峻程度。通过对此昆虫的基本了解，发觉气候因素对昆虫的

发生发展有着亲密关系，可以干脆影响昆虫的生长、发育、生存、繁殖,

从而造成害虫不同的发生期、发生量和危害程度。同时水稻螟虫是变温昆

虫，其生长、发育和繁殖与气象条件的关系极为亲密，所以我们从气候因

素角度入手进行分析是合理的。

本文从影响害虫生存繁殖的气候因素角度入手，结合往年的气象资料

以与影响害虫生存繁殖的重要气候因素，选取了平均气温、最低气温、日

照时间与降雨量四个主要影响因素，运用统计学方法确定虫害的发生量与

气候因子的关系，并给出相应的预料方法。

三.方法简介

针对实际问题，我们须要通过对历史数据的分析，给出准则：当给定

新时期下每一样本对应的各项气候指标时，能精确的推断其对应的虫害程

度。下面介绍两种方法。

（一）.判别分析

判别分析是在分类确定的条件下，依据某一探讨对象的各种特征值判

别其类型归属问题的一种多变量统计分析方法，是用以判别个体所属群体

的一种统计方法。依据判别中的组数，可以分为两组判别分析和多组判别

分析；依据判别函数的形式，可以分为线性判别和非线性判别；依据判别

式处理变量的方法不同，可以分为逐步判别、序贯判别；依据判别标准不

同，可以分为距离判别、Fisher判别、Bayes判别法。

本文选用的是Fisher判别法。Fisher判别，是依据线性Fisher函

数值进行判别，运用此准则要求各组变量的均值有显著性差异。该方法的

基本思想是投影，即将原来在R维空间的自变量组合投影到维度较低的D

维空间去，然后在D维空间中再进行分类。投影的原则是使得每一类的差

异尽可能小，而不同类间投影的离差尽可能大。Fisher判别的优势在于对

分布、方差等都没有任何限制，应用范围比较广。另外，用该判别方法建

立的判别方差可以干脆用手工计算的方法进行新样品的判别，这在很多时

候是特别便利的。

（二）.BP神经网络

BP神经网络是1986年由Rumelhart和McCelland为首的科学家

小组提出，是一种按误差逆传播算法训练的多层前馈网络，是目前应用最

广泛的神经网络模型之一。BP网络能学习和存贮大量的输入-输出模式映

射关系，而无需事前揭示描述这种映射关系的数学方程。它的学习规则是

运用梯度下降法，通过反向传播来不断调整网络的权值和阈值，使网络的

误差平方和最小。BP神经网络模型拓扑结构包括输入层、隐藏层和输出

层。

BP神经网络，由信息的正向传播和误差的反向传播两个过程组成。

输入层各神经元负责接收来自外界的输入信息，并传递给中间层各神经元;

中间层是内部信息处理层，负责信息变换，依据信息变更实力的需求，中

间层可以设计为单隐层或者多隐层结构；最终一个隐层传递到输出层各神

经元的信息，经进一步处理后，完成一次学习的正向传播处理过程，由输

出层向外界输出信息处理结果。当实际输出与期望输出不符时，进入误差

的反向传播阶段。误差通过输出层，按误差梯度F降的方式修正各层权值，

向隐层、输入层逐层反传。周而复始的信息正向传播和误差反向传播过程，

是各层权值不断调整的过程，也是神经网络学习训练的过程，此过程始终

进行到网络输出的误差削减到可以接受的程度，或者预先设定的学习次数

为止。

本文运用的是单层前馈网络模型，一般称为三层前馈网或三层感知器,

即：输入层、中间层（也称隐层）和输出层。如下图所示它的特点是：各

层神经元仅与相邻层神经元之间相互全连接，同层内神经元之间无连接,

各层神经元之间无反馈连接，构成具有层次结构的前馈型神经网络系统。

四.数据处理与分析

数据来自于1980〜1988年6〜10月安徽庐江的月平均温度、月最

低温度、月日照时间、月降雨量以与水稻的月平均虫害程度，且数据都是

经过标准化之后的。

（一）.数据预处理

首先，将降雨量的最终一个数据中的空格消退，再把原始数据的前两

列合并为一列，并保存为CSV（逗号分隔）格式，之后用函数将数据读入

到R的内存中。

>shuju<-read.csv（'chongqing.csv,）

为了便利之后的操作，对各列进行命名，其中date代表年度时间，

xl〜x4分别代表平均气温、最低气温、日照时间、降雨量，y代表虫害程

度。同时给出各类别下的样本量个数。由于本文不接受时间序列预料的方

法进行分析，故第一列的数据运用不到，最终确定运用的数据集为剔除原

数据框第一列之后的数据框。从原始数据中可知样本量为54,不是很多，

假如再进行分块，可能由于样本量的限制无法得到最好的结果，因此，下

文中的训练集和测试集均选用原始数据来充当。

>names(shuju)<-c('date,,'xr,'x2','x31,,x4,,'y,)

>x<-shuju[,2:6]

>table(y)

y

1234

356211

将数据框X中的变量链接到内存中，降低后续代码的困难程度，提高

代码的可读性。

>attach(x)

(二).判别分析

1.模型建立

安装并加载软件包MASSo

>install.packages(^MASS^)

>library(MASS)

将训练集中的y变量作为判别变量，其他剩余的变量作为特征变量,

运用lda()函数来进行费希尔判别下的线性判别。

>lda<-lda(y-.,data=x)#线性判别

>Ida

Call:

lda(y〜data=x)

Priorprobabilitiesofgroups:

1234

0.648148150.111111110.037037040.20370370

Groupmeans:

xlx2x3x4

1-0.03336857-0.007231429-0.0008800-0.2286971

20.769233330.7488000000.1302167-0.2696333

30.377650000.4366000000.6406500-0.4973500

40.565145450.6081909090.2076455-0.2936727

Coefficientsoflineardiscriminants:

LD1LD2LD3

xl1.1966105-1.8868126-3.1952390

x20.97735881.18404273.3471367

x30.24311362.1779657-0.2219837

x4-0.2754396-0.37820630.4283132

Proportionoftrace:

LD1LD2LD3

0.87670.11970.0036

输出结果的其次项为此次过程中各类别所运用的先验概率，其实就是

训练集中判别变量各水平下的数量占总样本量的比例。第三项为各变量在

每一类别中的均值，从中可以看出X1〜x3在各类别下有较明显的差别，

相反x4在各类别下的差别很小，也就说明降雨量对虫害程度的影响不是

很显著。第四项给出线性判别式的参数矩阵。第五项给出各线性判别式分

别的贡献比例。

>plot(lda)

从上图可以看到，在3个线性判别式下1〜4这4个类别的分布状况，

不同类别样本已经用相应数字标出。可能受样本量的限制，图中并不能看

出很明显的分布状况。

下面用上述模型对测试集进行预料，并绐出测试集中y变量的预料结

果与实际类别的混淆矩阵。

>pred<-predict(Ida,x[,k4])

>pred$class

[1]1144111144111111111114111

1441111441111421

[42]1114411111111

Levels:1234

>table(y,pred$class)

y1234

132003

22103

31001

46005

上述混淆矩阵的行表示实际的类别，列表示预料判定的类别。在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(x)

>e」da

[1]0.2962963

通过计算，上述模型的误判率为29.6%,可以看出该模型的预料效果

并不好，因此须要进一步的改进。

2.模型优化

从上文中我们知道，运用Fisher判别要求各组变量的均值有显著性

差异，从上述结果已经知道降雨量在各类别下的均值差异性很小，所以我

们可以考虑将其剔除，通过计算各变量在各类别下均值的方差，来作为剔

除某一变量的理论依据。

>var=>c(var(lda$means[,l]),var(lda$means[,2]),

+var(lda$means[,3]),var(lda$means[,4]))

>var

[1]0.116789090.107833140.077188360.01433243

从结果来看，x4在各类别下均值的方差相对其他三者而言差别较大,

因此考虑将其剔除，然后重新进行线性判别分析。

>lda2<-Ida(y-x1+x2+x3,data=x)#线性判别

>pred2<-predict(lda2,x[,1:3])

>table(y,pred2$class)

y1234

133002

23003

31001

45006

>e_lda2<-sum(as.numeric(pred2$class)!=y)/nrow(x)

>e_lda2

[1]0,2777778

从误判率的角度来看，剔除x4之后的模型误判率为27.8%,相对原

来的模型有所降低；但从混淆矩阵来看，对于第2类和第3类的误判率为

100%,这相对原来的模型更难让人接受。两个模型各有好坏，但整体而

言，都不是很好。

3.模型应用

假如运用上述两个模型对将来进行预料的话，只需将所须要的变量存

放至一个新的数据框，须要留意的是，数据框中各变量的名字须要和建立

模型时运用的数据框中各变量的名称一样，即平均气温、最低气温、日照

时间、降雨量分别用xl〜x4来表示。然后干脆运用模型来给出预料结果。

(三).BP神经网络

1.模型建立

安装并加载软件包nncto

>install.packages(“nnet")

>library(nnet)

首先在建模之前，先介绍软件包nnet中的一个函数：class.ind()o

该函数可以通过类别变量的因子变量来生成一个类指标矩阵。下面通过一

个简洁的例子来更清晰地表明该函数的功能。

>z=c(l,2,3)#定义一个水平数为3的向量

>class.ind(z)

123

[1,]100

[2,]010

[3,]001

从输出结果可以看到，该函数主要是将向量变成一个矩阵，其中每行

代表一个样本。只是将样本的类别用。和1来表示，即假如是该类，则在

该类别名下用1表示，而其余的类别名下面用0表示。

将数据框的前4列作为自变量数据,class.ind(y)作为类别变量数据，

设置隐藏层节点个数为4,模型的最大迭代次数为10。。次，用来防止模

型的死循环。之后用核心函数nnet()来建立神经网络模型。

>model=nnet(x[1:4],class.ind(y),data=x,size=4,

+decay=5e-4,maxit=1000)#神经网络

>P=x[,l: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=max.col(pred3)#确定测试集中每一个样本所属类别

>table(y,pnew)#给出混淆矩阵

pnew

y1234

133101

22400

30011

410010

>e_e<-sum(pnew!=y)/nrow(x)#计算误判率

>e_e

[1]0.1111111

从结果可以看出，该模型的误判率为11.1%,相比较判别分析大幅度

降低，说明该模型对测试集的预料效果很好。

>names(model)#模型中包含的全部输出结果

[1]"n""nunits""nconn""conn"nnsunits""decay"[7]

"entropy"usoftmax""censored""value""wts"

"convergence"

[13]"fitted.values""residuals""call"

其中：value表示迭代最终值。

wts表示模型在迭代过程中所找寻到的最优权重值，其取值为系统随

机生成，即每次建模所运用的迭代初始值都是不相同的。由此可以看出神

经网络模型的敏捷性很高。

convergence表示在模型建立的迭代过程中，迭代次数是否达到最大

迭代次数。假如结果为1,则表明达到最大迭代次数；假如结果为。，则

表明没有达到最大迭代次数。假如结果没有达到最大迭代次数，说明该模

型可以被接受，否则，就要增大参数maxit的设定值，重新建模。

>model$convergence

[1]0

该模型的convergence取值为0,说明该模型可以被接受。

2.模型优化

虽然上述模型的各项指标可能都满意要求，但是，隐藏层中节点的个

数和迭代次数都是人为设定的，随着两者数值的变更建立的模型也不尽相

同，所以我们须要找寻两者的最优值，使其对应的模型最优。下面通过编

写函数来实现这一目的。

(1).找寻隐藏层节点个数最优值

>f=function(n){

+e=matrix(02n,dimnames=list(c('误判概率

+'迭代终值工c(l:n)))#定义一矩阵，用于存放结果

+for(iinl:n){

+model=nnet(x[1: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)

+e[l,i]=e_e

+e[2,i]=model$value

+1

+print(e)

+print(plot(as.numeric(e[l,]),type='r,col='red',

+main=候判率变更曲线3ab±隐藏层结点个数;ylab

+='误判率’))

+}

>f(10)#将。取为10

输出结果:

1234

5误判概率0.25925930.20370370.18518520.07407410.0

370371

迭代终值20.70701215.36959815.7799998.23342796.50628

53

67891

0

误判概率0.01851850.03703710.0000000.03703710.01851

85

迭代终值4.49751954.66498322.7332964.39859893.29612

55

误判率变化曲线

树卜

声:

呸

结合矩阵以与折线图，可以看出当n=8的时候，误判率最小，同时迭

代最终值也达到最小，说明当n=8即隐藏层节点个数为8时，模型的拟

合效果最好，对测试集做出的预料效果也最好。

⑵.找寻迭代次数最优值

>ff=function(n.l{

+e=0

+for(iinl:n){

+model=nnet(x[,1-4],class.ind(x$y),data=x,

+size=8,decay=5e-4,maxit=i)

+e[i]=sum(max.col(predict(model,x[,1:4]))!=

+x$y)/nrow(x)

+}

+print(plot(l:n,e,T,ylab=，误判概率',xlab='训练周

+期',col='blue',ylim=c(min(e),max(e))))

ff(1000)

训练周期

从模型对测试集的误判率随训练周期的变更曲线可以看出，误判概率

随着训练周期的增大先是急速下降，后趋于平稳。但是也有出现反向增长

的现象，可能与训练集和测试集的选取有关。不过，从图中整体的变更趋

势可以看出曲线的拐点在maxit=200处，所以maxit的最优值可能为2

OOo

3.建模重建

取size=8,maxit=200,重新建立神经网络模型，得到如下结果：

>modell=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_el<-sum(pnewl!=y)/nrow(x)

>e_el

[1]0,05555556

>model1$convergence

[1]1

误判概率为5.6%,较上次模型有所下降，但是模型的迭代次数达到

了最大迭代次数，因此该模型可能还不是最优的。

取size=8,maxit=1000,重新建立神经网络模型，得到如下结果:

>model2=nnet(x[1:4],class.ind(y),data=x,size=8,

+decay=5e-4,maxit=1000)

>prcd5=prcdict(modcl2,p)

>pnew2=max.col(pred5)

>table(y,pnew2)

Pnew2

y