施肥效分析MatlabWord文档下载推荐.docx
- 文档编号:19646914
- 上传时间:2023-01-08
- 格式:DOCX
- 页数:15
- 大小:51.97KB
施肥效分析MatlabWord文档下载推荐.docx
《施肥效分析MatlabWord文档下载推荐.docx》由会员分享,可在线阅读,更多相关《施肥效分析MatlabWord文档下载推荐.docx(15页珍藏版)》请在冰豆网上搜索。
参数估计值
参数置信区间
a0
-12.8361
[-20.6921,-4.9802]
a1
0.1903
[0.1597,0.2208]
a2
0.0842
[0.0418,01265]
a3
0.0735
[0.0512,0.0958]
a11
-0.0003
[-0.0004,-0.0003]
a22
-0.0002
[-0.0003,0]
a33
-0.0001
[-0.0001,0]
R^2=0.9190F=43.4925P<
<
0.0000S^2=6.1094
表1模型(4)的计算结果
结果分析表1表示,R2=0.9190指因变量y(产量)的91.90%可由模型确定,F值远远超过F检验的平均值,p值小于α,因此模型(4)整体来说是可用的。
从表1的回归系数给出模型(4)中a0,a1,a2,a3,a11,a22,a33的估计值,即a0=-12.8361,a1=0.1903,a2=0.0842,a3=0.0735,a11=-0.0003,a22=-0.0002,a33=-0.0001.检验它们的置信区间发现它们的置信区间都不包含零点,说明此模型中的回归变量对模型的影响都还是显著地,而此模型的建立也是基本合理的。
则此时
y=-12.8361+0.1903y1+0.0842y2+0.0735y3-0.0003y12-0.0002y22-0.0001y32
(5)
我们可以根据此模型来预测一定N,P,K施肥量下,所可能得到的产量的值。
3.2模型改进
模型(4)中回归变量y1,y2,y3对因变量y的影响是相互独立的,即土豆的产量y的均值与N,P,K各自的施肥量有关,而就我们对实际生活中的施肥问题,就生物,化学等机理上进行分析,其实很清楚元素间都会有相互作用的,有些元素共同作用就会出现相互促进或者相互抑制的作用,即我们可以猜想N,P,K它们之间的相互作用会对产量y有影响。
在此分析上模型可以改进为
y=a0+a1*y1+a2*y2+a3*y3+a11*y12+a22*y22+a33*y32+
a12*y1*y2+a13*y1*y3+a23*y2*y3+ε(6)
即在模型四的基础上增加了NP,NK,PK之间的相互影响。
同模型四的解法一样,我们用matlab求解模型(6)(具体程序见附件【2】)。
得到模型(6)的回归系数及其置信区间(置信水平α=0.05),检验统计量R2,F,p,s2的结果见表2.
参数估计值
[0,0]
[0,0]
0.1813
[0.1633,0.1992]
[-0.0004,-0.0003]
[-0.0003,0]
[-0.0001,0]
a12
0.001
[0.0009,0.0012]
a13
a23
-0.0005
[-0.0006,-0.0004]
R^2=0.9190F=43.4925P<
表2模型(6)的计算结果
表1和表2相比,发现参数的置信区间都不包含零点,但是回归模型的监测统计量都相同,特别的R2相同,可以发现模型(6)在模型(4)的基础上并没有得到更好地改进。
3.3进一步思考后改进模型
于是我们自然的联想到,当3种主要元素的化肥共同添加的到作物上时。
共同作用的可能并不仅仅是两两共同作用,而是3中元素都会相互作用。
于是我们想到添加含有3个因素变量的式子,即求当3中因素共同作用于作物时,对产量的影响。
此时模型可以表示为
y=a0+a1*y1+a2*y2+a3*y3+a11*y12+a22*y22+a33*y2+a12*y1y2+a13*y1y3
+a23*y2y3+a111*y13+a222*y23+a333*y3+a123*y1y2y3+ε(7)
即此模型在(6)的基础上除了还添加了NNN,PPP,KKK,NPK作用的影响。
其中NNN可以理解为N的量在作物上添加的很多。
同理可以理解PPP,KKK分别为P和K的量在作物上添加的很多。
而NPK可以理解为是氮,磷,钾共同作用于作物时对作物产量的影响。
另外,可以发现模型六回归分析中相关参数估计值很小,由于模型建立式中存在NNN相乘项,因此对于施肥量中的所有数据,我们进行一个处理。
即将每一组数据除以改组数据中最大值,使得施肥量不再是数值大小的关系而是一个比例关系。
这样可以减少对小参数的忽略问题,是问题更加精确得以解决。
同样用matlab求解的表3.
-15.3594
[-21.8224,-8.8965]
76.7818
[44.4492,109.1145]
38.7451
[6.4679,71.0224]
109.2842
[76.9844,141.5840]
-45.7793
[-123.7268,32.1681]
-49.2256
[-125.7797,27.3285]
-182.6483
[-259.2521,-106.0446]
[0,0]
a111
-14.9509
[-65.3347,35.4329]
a222
20.5324
[-28.4887,69.5535]
a333
99.7052
[50.6407,148.7697]
a123
R^2=0.9586F=51.4201P<
0.0000S^2=3.5934
表3模型(7)的计算结果
有表3可知,R2=0.9586指因变量y(产量)的95.86%可由模型确定,
这个值与表2中的0.9109相比有了更为准确的模型解。
F值远远超过F
检验的平均值,p值也远远小于α,并且s2也有模型(6)中的6.1094变
成3.5934,说明剩余方差更小。
因此模型(7)整体来说是可用的,并且
要明显优于模型(4)和模型(6)。
但模型(7)也有不足之处,我们发现回归系数中a11,a22包含零点,说明NN,PP对产量的影响很小。
而这也严重影响模型的准确率。
并且可以看到模型中a12,a13,a23,a123都为零。
说明这些相应变量几乎对模型没有影响。
这就是我们接下来改进的方向。
但对于接下来的分析,因为模型(7)要明显优于模型(4)和模型(6)。
我们我们重点在模型(7)的基础上进行进一步分析和改进。
图4模型(7)中农作物产量y与N,P,K之间的残插图
3.4、在3.3模型基础上去除异常数据后修正模型
因为回归系数中a11,a22包含零点,为了寻找改进的方法,我们采用残差分析的方法(残差ε值农场品产值的实际值y与用模型(7)估计得到的产值之差),得到模型(7)的农作物产量y与N,P,K的施肥量的残插图(具体matlab程序见附录【4】)。
图4就是模型(7)的农作物产量y与N,P,K的施肥量的残插图,其中红色的*代表N元素与y的关系,绿色的*代表P元素与y的关系,而蓝色的*代表K元素与y的关系,
从图四的残插图分析中,我们可以发现数据中有几个残差局势特别大的异常点,即当N肥施加量为101(公斤/公顷)时产量为32.29(吨/公顷);
p肥施加量为24(公斤/公顷)时产量为32.47(吨/公顷)和K肥施加量为279(公斤/公顷)时产量为37.73(吨/公顷)时,这三个点误差最大。
所以我们考虑把这3个点做适当的修改下,然后求出改进后的模型。
与元模型(7)做比较。
同样,我们看题目所给的N,P,K的数据,发现其实数据中存在很大的误差,如图6中红色数据所示,因为在单一变量的情况下,随着变量的增加,不可能出现自变量反复变化的情况。
同时我们应注意到,删除异常数据的过程中,我们要注意其删除后仍为一个矩阵,方便计算,同时也不能删除过多的数据,因此我们选择在N,P,K数据中各删除一组偏差比较大的数据。
N
P
K
施肥量
(公斤/公顷)
产量
(吨/公顷)
34
67
101
135
202
259
336
404
471
15.18
21.36
25.72
32.29
34.03
39.45
43.15
43.46
40.83
30.75
24
49
73
98
147
196
245
294
342
33.46
32.47
36.06
37.96
41.04
40.09
41.26
42.17
40.36
42.73
47
93
140
186
279
372
465
558
651
18.98
27.35
34.86
38.52
38.44
37.73
38.43
43.87
42.77
46.22
图6题目数据中的异常反应
我们考虑到这可能是做此调查实验时或许出现太阳暴晒导致化肥挥发或者持续雨水天气导致化肥流失的情形,从而导致当时采样数据的不准确,而影响模型(7)的准确性。
正如对比图4中模型(7)的残差图一样,同样可以发现在题目数据异常的地方,残差都是比较大的。
这些都给我们修改异常数据提供了理论和实际的支持。
-15.44866
[-21.0653,-9.8319]
80.7213
[52.2631,109.1795]
31.1933
[3.6104,58.7763]
115.3799
[86.5207,144.2391]
-51.1202
[-116.0380,13.7976]
-36.3340
[-97.1984,24.5305]
-192.9619
[-256.7614,-129.1263]
-13.5669
[-54.5730,27.4393]
14.1670
[-24.2120,52.5460]
104.2762
[64.6133,143.9391]
R^2=0.9730F=67.9596P<
0.0000S^2=2.7107
图67删除异常数据后修正模型回归分析图
而在修改数据的基础上,我们发现R2=0.9730指因变量y(产量)的967.30%可由模型确定,这个值与表3中的0.9586相比有了更为准确的模型解。
F值远远超过F检验的平均值,p值也远远小于α,并且s2也有模型(6)中的3.5934变成2.7107,说明剩余方差更小。
因此,可以说明修改数据后模型的准确率得到的提高,从而新模型要优于模型(7)。
图8删除异常数据后的修正模型残差图
对于修正模型后的残差图我们也可以看到残差值平均值的期望接近于零,比修改前的模型更为合理,因此我们认为这时建立的模型比较可靠。
4、附录:
附录【1】模型(4)解的matlab程序
%y=a0+a1*N+a2*P+a3*K+a11*NN+a22*PP+a33*KK
d1=[0.000034.000067.0000101.0000135.0000202.0000259.0000336.0000404.0000471.0000];
d2=[0.000024.000049.000073.000098.0000147.0000196.0000245.0000294.0000342.0000];
d3=[0.000047.000093.0000140.0000186.0000279.0000372.0000465.0000558.0000651.0000];
y=[15.180021.360025.720032.290034.030039.450043.150043.460040.830030.7500...
33.460032.470036.060037.960041.040040.090041.260042.170040.360042.7300...
18.980027.350034.860038.520038.440037.730038.430043.870042.770046.2200]'
;
d1=[d1,d1(7)*ones(1,10),d1(7)*ones(1,10)]'
d2=[d2(7)*ones(1,10),d2,d2(7)*ones(1,10)]'
d3=[d3(7)*ones(1,10),d3(7)*ones(1,10),d3]'
x=[ones(30,1),d1,d2,d3,d1.*d1,d2.*d2,d3.*d3];
[b,bint,r,rint,stats]=regress(y,x);
附录【2】模型(6)解的matlab程序
%y=a0+a1*N+a2*P+a3*K+a11*NN+a22*PP+a33*KK+a12*NP+a13*NK+a23*PK
x=[ones(30,1),d1,d2,d3,d1.*d1,d2.*d2,d3.*d3,d1.*d2,d1.*d3,d2.*d3];
附录【3】模型(7)解的matlab程序
%y=a0+a1*N+a2*P+a3*K+a11*NN+a22*PP+a33*KK+a12*NP+a13*NK+a23*PK+a111*NNN+a222*PPP+a333*KKK+a123*NPK;
d1=[03467101135202259336404471]./471;
d2=[024497398147196245294342]./342;
d3=[04793140186279372465558651]./651;
y1=[15.180021.360025.720032.290034.030039.450043.150043.460040.830030.7500];
y2=[33.460032.470036.060037.960041.040040.090041.260042.170040.360042.7300];
y3=[18.980027.350034.860038.520038.440037.730038.430043.870042.770046.2200];
y=[y1y2y3]'
x=[ones(30,1),d1,d2,d3,d1.*d1,d2.*d2,d3.*d3,d1.*d2,d1.*d3,d2.*d3,d1.*d1.*d1,d2.*d2.*d2,d3.*d3.*d3,d1.*d2.*d3];
附录【4】模型(7)的相关残插图
%y=a0+a1*N+a2*P+a3*K+a11*NN+a22*PP+a33*KK+a12*NP+a13*NK+a23*PK+
a111*NNN+a222*PPP+a333*KKK+a123*NPK;
symsxyzw
y=0.5731;
z=0.5714;
x=[03467101135202259336404471]./471;
%w1=0.1813*z-0.0003.*x.^2-0.0002.*y.^2-0.0001*z^2+0.001.*x.*y-0.0005.*y*z;
w1=76.7818.*x+38.7451.*y+109.2842.*z-45.7793.*x.^2-49.2256.*y.^2-182.6483.*z.^2-14.9509.*x.^3+20.5324.*y.^3+99.7052.*z.^3-15.3594;
x=0.5499;
z=0.5714;
y=[024497398147196245294342]./342;
%w2=0.1813*z-0.0003.*x.^2-0.0002.*y.^2-0.0001*z^2+0.001.*x.*y-0.0005.*y*z;
w2=76.7818.*x+38.7451.*y+109.2842.*z-45.7793.*x.^2-49.2256.*y.^2-182.6483.*z.^2-14.9509.*x.^3+20.5324.*y.^3+99.7052.*z.^3-15.3594;
y=0.5731;
z=[04793140186279372465558651]./651;
%w3=0.1813.*z-0.0003.*x.^2-0.0002.*y.^2-0.0001.*z.^2+0.001.*x.*y-0.0005.*y.*z;
w3=76.7818.*x+38.7451.*y+109.2842.*z-45.7793.*x.^2-49.2256.*y.^2-182.6483.*z.^2-14.9509.*x.^3+20.5324.*y.^3+99.7052.*z.^3-15.3594;
%氮肥对产量影响残差图
ch1=[15.180021.360025.720032.290034.030039.450043.150043.460040.830030.7500];
plot(ch1,w1-ch1,'
r*'
)
holdon
%磷肥对产量影响残差图
ch2=[33.460032.470036.060037.960041.040040.090041.260042.170040.360042.7300];
plot(ch2,w2-ch2,'
g*'
%钾肥对产量影响残差图
ch
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 施肥 分析 Matlab
![提示](https://static.bdocx.com/images/bang_tan.gif)