2.4 建立模型
Yp表示施磷肥时土豆的产量，tp表示磷肥的施肥量，模型设为
【3.2】
取Y=ln(44-Yp),a0=lnB,a1=-C,【3.2】就等价于下面的线性回归模型
【3.3】
2.5 参数求解
matlab计算代码如下
>> clear >> load d:\tudou >> yp(2)=(yp(1)+yp(3))/2; >> tp1=(tp-mean(tp))/std(tp); >> Y=log(44-yp); >> plot(tp1,Y,'*') >> [mean(tp),std(tp)] ans = 146.8000 118.2547 >> [a,at,r,rt,st]=regress(Y',[ones(length(tp),1),tp1']); >> a a = 1.4041 -0.6211 >> at at = 1.1495 1.6587 -0.8895 -0.3527 >> st st = 0.7807 28.4829 0.0007 0.1219 >> rt rt = -0.5561 0.9161 -0.5894 0.9383 -0.6351 0.9434 -0.8102 0.8236 -1.2249 0.0744 -0.8761 0.7971 -0.9568 0.6814 -1.0436 0.4756 0.1942 1.1279 -0.8104 0.5307
第9个样本异常，修正后计算
>> Y(9)=[];tp1(9)=[]; [a,at,r,rt,st]=regress(Y',[ones(length(tp1),1),tp1']); >> st st = 0.9154 75.7884 0.0001 0.0536 （极显著） a = 1.3133 -0.7467
回归结果为
实验值与理论值对比如下图所示。
理论值 实验值
34.6242 33.4600
35.9398 34.7600
37.1144 36.0600
38.0806 37.9600
38.9432 41.0400
40.2863 40.0900
41.2726 41.2600
41.9970 42.1700
42.5290 40.3600
42.9129 42.7300
3、土豆产量与磷肥关系（2）
在上面的土豆产量与磷肥关系（1）里，我们假定土豆的最大产量为44后才得以线性化（这个假设有时候不太合理），如果不能线性化，就只能非线性拟合。
3.1 经验模型
3.2 模型求解
调用matlab的最小二乘曲线拟合函数lsqcurvefit,来求解模型参数beta。先编写模型的经验函数的m文件tpfun.m。
调用格式为：
>> clear >> load d:\tudou >> x=tp';y=yp'; >> [beta,res]=lsqcurvefit(@tpfun,[40,0.2,2],x,y)
[beta,res]=lsqcurvefit(@tpfun,[40,0.2,2],x,y)
- beta为返回的参数的估计值，按照tpfun.m中顺序给出
- res 返回残差平方和，即经验值与实验值之间的误差平方和
- @tpfun 调用函数名为tpfun.m的m文件
- [40,0.2,2] 是参数beta的初始值
- X 土豆试验中的磷肥施肥量（列向量）
- Y 土豆实验中的土豆产量(列向量）
计算结果为
beta = 19.4643 3.9464 27.4523 res = 18.1411
即土豆对磷肥的拟合方程为：
3.3 误差分析
土豆对磷肥的拟合误差为
经验值与实验值对比图像
实验值y 回归值y* 33.4600 32.5365 32.4700 35.0157 36.0600 36.5785 37.9600 37.6559 41.0400 38.5330 40.0900 39.8342 41.2600 40.8111 42.1700 41.5935 40.3600 42.2462
土豆产量的实验值与拟合值差的分布图
>> E=y-tpfun(beta,x) E = 0.9235 -2.5457 -0.5185 0.3041 2.5070 0.2558 0.4489 0.5765 -1.8862 -0.0654
3.4 正态性假设检验
（1）绘制正态概率图
>> normplot(E)
除了3个点外，其余的点都均匀分布在直线两侧，即数据近似服从正态分布
(2)Jbtest正态检验
>> [h,p,jbstat,critval]=jbtest(E,0.05) h = 0 p = 0.5000 jbstat = 0.1221 critval = 2.5239
Jbtest检验：h=1时拒绝原假设H0；当jbstat>critval时拒绝原假设。计算结果显示，没有理由拒绝H0，即接受H0。
认为E服从均值为0的正态分布。
