施肥效分析Matlab.doc

施肥效果分析模型一、摘要（略）二、问题的分析为了大致的了解氮(N)、磷(P)、钾(K)对施肥效果的影响，首先我们利用表中的数据用matlab分别做出了产量对氮(N)、磷(P)、钾(K)的散点图（见图1，图2和图3中的*号）。 图1,N与产量的关系 图2，K与产量的关系 图3，P和产量的关系从3个图中我们可以发现氮(N)、磷(P)、钾(K)对产量的影响都不是简单线性关系，从图1中我们可以明显的看到N对产量的影响是先增加后减少，呈现二次函数的关系，于是我们很自然用二次函数的模型来表示N的影响。至于图2和图3中P，K对作物产量的影响上，可能由于数据范围小和数据过少等原因，我们并不能明确的从图中推算P,K对产量是呈现怎么的趋势，但是可以明确的是它们的关系绝不是简单的线性关系，我们姑且可以把它当做二次函数的模型来求解，并且我们根据图2，图3的数据分布，可以推测这其中有些数据存在很大的误差，在以后的模型分析中可以考虑将这些异常点去掉。3、 模型的建立 3.1基本模型 记土豆的产量是y（吨/公顷），N的使用量为y1（公斤/公顷），P的使用量为y2（公斤/公顷），K的使用量为y3（公斤/公顷），基于上面的分析，我们可以建立初步模型。N对产量的影响可以表示为 y=a0+a1*y1+a11*y12 + （1）同理，P,K对产量的影响可以表示为 y=a0+a1*y2+a22*y22 + （2）综合以上分析，结合模型（1）（2）和（3）可以建立一下模型。 可得 y=a0+a1*y1+a2*y2+a3*y3+a11*y12+a22*y22+a33*y32 + （4） 其中（4）式右端的y1，y2，y3称为回归变量。而a0，a1，a2，a3，a11，a22，a33称为回归系数。由题目表中的数据我们估计，影响y的其他因素都包含在随即误差中，如果模型选择的合适，应大致符合均值为0的正太分布。这个在之后的证明中会检验。 3.1.2模型求解 对于此模型，我们可以直接用matlab统计工具箱的命令regress来求解（具体程序见附录【1】）。其中输出b为回归系数的估计值，而bint为b的置信区间，stats为回归模型的检验统计量，有四个值，第1个是回归方程的决定系数R2，第2个是F统计量，第3个是F统计量对应的概率值p，第4个是剩余方差s2。得到模型（4）的回归系数及其置信区间（置信水平=0.05），检验统计量R2，F，p，s2 的结果见表1. 参数 参数估计值 参数置信区间 a0-12.8361-20.6921,-4.9802 a10.19030.1597,0.2208 a20.08420.0418,01265 a30.07350.0512,0.0958 a11-0.0003-0.0004,-0.0003 a22-0.0002-0.0003,0 a33-0.0001-0.0001,0R2=0.9190 F=43.4925 P0.0000 S2=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. 参数参数估计值 参数置信区间 a00 0 ,0 a10 0 ,0 a20 0 ,0 a30.1813 0.1633,0.1992 a11-0.0003 -0.0004,-0.0003 a22-0.0002 -0.0003,0 a33-0.0001 -0.0001,0 a120.001 0.0009,0.0012 a130 -0.0001,0 a23-0.0005 -0.0006,-0.0004 R2=0.9190 F=43.4925 P0.0000 S2=6.1094 表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. 参数参数估计值 参数置信区间 a0 -15.3594 -21.8224,-8.8965 a176.7818 44.4492 ,109.1145 a238.7451 6.4679 ,71.0224 a3109.2842 76.9844,141.5840 a11-45.7793 -123.7268,32.1681 a22-49.2256 -125.7797,27.3285 a33-182.6483 -259.2521,-106.0446 a120 0,0 a130 0,0 a230 0,0 a111-14.9509 -65.3347,35.4329 a22220.5324 -28.4887,69.5535 a33399.7052 50.6407,148.7697 a1230 0,0 R2=0.9586 F=51.4201 P0.0000 S2=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数据中各删除一组偏差比较大的数据。NPK施肥量（公斤/公顷）产量（吨/公顷）施肥量（公斤/公顷）产量（吨/公顷）施肥量（公斤/公顷）产量（吨/公顷）0346710113520225933640447115.1821.3625.7232.2934.0339.4543.1543.4640.8330.7502449739814719624529434233.4632.4736.0637.9641.0440.0941.2642.1740.3642.730479314018627937246555865118.9827.3534.8638.5238.4437.7338.4343.8742.7746.22 图6 题目数据中的异常反应 我们考虑到这可能是做此调查实验时或许出现太阳暴晒导致化肥挥发或者持续雨水天气导致化肥流失的情形，从而导致当时采样数据的不准确，而影响模型（7）的准确性。正如对比图4中模型（7）的残差图一样，同样可以发现在题目数据异常的地方，残差都是比较大的。这些都给我们修改异常数据提供了理论和实际的支持。 参数参数估计值 参数置信区间 a0 -15.44866 -21.0653,-9.8319 a180.7213 52.2631 ,109.1795 a231.1933 3.6104 ,58.7763 a3115.3799 86.5207,144.2391 a11-51.1202 -116.0380,13.7976 a22-36.3340 -97.1984,24.5305 a33-192.9619 -256.7614,-129.1263 a120 0,0 a130 0,0 a230 0,0 a111-13.5669 -54.5730,27.4393 a22214.1670 -24.2120,52.5460 a333104.2762 64.6133,143.9391 a1230 0,0 R2=0.9730 F=67.9596 P0.0000 S2=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*KKd1=0.0000 34.0000 67.0000 101.0000 135.0000 202.0000 259.0000 336.0000 404.0000 471.0000;d2=0.0000 24.0000 49.0000 73.0000 98.0000 147.0000 196.0000 245.0000 294.0000 342.0000;d3=0.0000 47.0000 93.0000 140.0000 186.0000 279.0000 372.0000 465.0000 558.0000 651.0000;y=15.1800 21.3600 25.7200 32.2900 34.0300 39.4500 43.1500 43.4600 40.8300 30.7500 . 33.4600 32.4700 36.0600 37.9600 41.0400 40.0900 41.2600 42.1700 40.3600 42.7300 . 18.9800 27.3500 34.8600 38.5200 38.4400 37.7300 38.4300 43.8700 42.7700 46.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*PKd1=0.0000 34.0000 67.0000 101.0000 135.0000 202.0000 259.0000 336.0000 404.0000 471.0000;d2=0.0000 24.0000 49.0000 73.0000 98.0000 147.0000 196.0000 245.0000 294.0000 342.0000;d3=0.0000 47.0000 93.0000 140.0000 186.0000 279.0000 372.0000 465.0000 558.0000 651.0000;y=15.1800 21.3600 25.7200 32.2900 34.0300 39.4500 43.1500 43.4600 40.8300 30.7500 . 33.4600 32.4700 36.0600 37.9600 41.0400 40.0900 41.2600 42.1700 40.3600 42.7300 . 18.9800 27.3500 34.8600 38.5200 38.4400 37.7300 38.4300 43.8700 42.7700 46.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,d1.*d2,d1.*d3,d2.*d3;b,bint,r,rint,stats=regress(y,x);附录【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=0 34 67 101 135 202 259 336 404 471./471;d2=0 24 49 73 98 147 196 245 294 342./342;d3=0 47 93 140 186 279 372 465 558 651./651;y1=15.1800 21.3600 25.7200 32.2900 34.0300 39.4500 43.1500 43.4600 40.8300 30.7500;y2=33.4600 32.4700 36.0600 37.9600 41.0400 40.0900 41.2600 42.1700 40.3600 42.7300 ;y3=18.9800 27.3500 34.8600 38.5200 38.4400 37.7300 38.4300 43.8700 42.7700 46.2200;y=y1 y2 y3;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,d1.*d2,d1.*d3,d2.*d3, d1.*d1.*d1,d2.*d2.*d2,d3.*d3.*d3,d1.*d2.*d3;b,bint,r,rint,stats=regress(y,x);附录【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;syms x y z wy=0.5731;z=0.5714;x=0 34 67 101 135 202 259 336 404 471./471;%w1=0.1813*z-0.0003.*x.2-0.0002.*y.2-0.0001*z2+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=0 24 49 73 98 147 196 245 294 342./342;%w2=0.1813*z-0.0003.*x.2-0.0002.*y.2-0.0001*z2+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;x=0.5499; y=0.5731;z=0 47 93 140 186 279 372 465 558 651./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.1800 21.3600 25.7200 32.2900 34.0300 39.4500 43.1500 43.4600 40.8300 30.7500;plot(ch1,w1-ch1,r*) hold on%磷肥对产量影响残差图ch2=33.4600 32.4700 36.0600 37.9600 41.0400 40.0900 41.2600 42.1700 40.3600 42.7300;plot(ch2,w2-ch2,g*) hold on%钾肥对产量影响残差图ch3=18.9800 27.3500 34.8600 38.5200 38.4400 37.7300 38.4300 43.8700 42.7700 46.2200;plot(ch3,w3-ch3,b*)legend(N,P,K,Location,northwest);附录【5】模型（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=0 34 67 135 202 259 336 404 471./471;d2=0 49 73 98 147 196 245 294 342./342;d3=0 47 93 140 186 372 465 558 651./651;y1=15.1800 21.3600 25.7200 34.0300 39.4500 43.1500 43.4600 40.8300 30.7500;y2=33.4600 36.0600 37.9600 41.0400 40.0900 41.2600 42.1700 40.3600 42.7300 ;y3=18.9800 27.3500 34.8600 38.5200 38.4400 38.4300 43.8700 42.7700 46.2200;y=y1 y2 y3;d1=d1,d1(7)*ones(1,9),d1(7)*ones(1,9);d2=d2(7)*ones(1,9),d2,d2(7)*ones(1,9);d3=d3(7)*ones(1,9),d3(7)*ones(1,9),d3;x=ones(27,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;b,bint,r,rint,stats=regress(y,x);残差分析图源程序% 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;syms x y z wy=0.5731;z=0.5714;x=0 34 67 135 202 259 33