nlinfit+regress_data fitting.rar_作物模型_农作物_数据拟合regress+nlinfit

clc clear ETm1=148.1;ETm2=111.8;ETm3=124.7;ETm4=89.4;Ym=7138.5; ETa1=[148.1 113.2 107.6 133.9 132.1 128.2 129.7 140.5 135.3 110.6 128.4 130.1]; ETa2=[111.8 96.6 88.3 91.0 77.9 99.4 92.5 112.9 108.0 83.3 90.4 102.6]; ETa3=[124.7 92.1 84.9 106.9 93.9 85.3 71.9 108.6 101.7 95.2 83.4 94.7]; ETa4=[89.4 67.2 64.2 70.3 65.0 78.7 69.4 68.6 65.0 72.3 73.6 61.4]; Ya=[7138.5 5757.0 4576.5 6111.0 4555.5 5520.0 5329.5 6345.0 6040.5 5076.0 5442.0 6130.5]; ET1=ETa1/ETm1; ET2=ETa2/ETm2; ET3=ETa3/ETm3; ET4=ETa4/ETm4; %% 加法模型1-Blank模型 X11=ET1;X12=ET2;X13=ET3;X14=ET4; Y11=Ya/Ym; T1=[X11' X12' X13' X14']; [b1,bint1,r1,rint1,stats1]=regress(Y11',T1) figure; subplot(3,2,1) plot(Ya,'k-') % Y12=Ym.*(0.0525*(ET1)+0.5575*(ET2)+0.2909*(ET3)+0.0607*(ET4)); % hold on;plot(Y12,'b-') Y13=Ym.*(b1(1)*(ET1)+b1(2)*(ET2)+b1(3)*(ET3)+b1(4)*(ET4)); hold on;plot(Y13,'r-') title('Blank模型') xlabel('蒸发蒸腾量(mm)') ylabel('产量(kg/hm^2)') % legend('真实值','函数拟合值') %% 加法模型2-Stewart模型 X21=(ETm1-ETa1)/ETm1;X22=(ETm2-ETa2)/ETm2;X23=(ETm3-ETa3)/ETm3;X24=(ETm4-ETa4)/ETm4; Y21=1-Ya/Ym; T2=[X21' X22' X23' X24']; [b2,bint2,r2,rint2,stats2]=regress(Y21',T2) % figure; subplot(3,2,3) plot(Ya,'k-') Y22=Ym.*(1-(b2(1)*((ETm1-ETa1)/ETm1)+b2(2)*((ETm2-ETa2)/ETm2)+b2(3)*((ETm3-ETa3)/ETm3)+b2(4)*((ETm4-ETa4)/ETm4))); hold on;plot(Y22,'r-') title('Stewart模型') xlabel('蒸发蒸腾量(mm)') ylabel('产量(kg/hm^2)') % legend('真实值','函数拟合值') %% 加法模型3-Singh模型 X31=1-(1-ET1).^2;X32=1-(1-ET2).^2;X33=1-(1-ET3).^2;X34=1-(1-ET4).^2; Y31=Ya/Ym; T3=[X31' X32' X33' X34']; [b3,bint3,r3,rint3,stats3]=regress(Y31',T3) % figure; subplot(3,2,5) plot(Ya,'k-') Y32=Ym.*(b3(1)*(1-(1-ET1).^2)+b3(2)*(1-(1-ET2).^2)+b3(3)*(1-(1-ET3).^2)+b3(4)*(1-(1-ET4).^2)); hold on;plot(Y32,'r-') title('Singh模型') xlabel('蒸发蒸腾量(mm)') ylabel('产量(kg/hm^2)') % legend('真实值','函数拟合值') %% 乘法模型1-Jensen模型 X41=log(ET1);X42=log(ET2);X43=log(ET3);X44=log(ET4); Y41=log(Ya/Ym); T4=[X41' X42' X43' X44']; [b4,bint4,r4,rint4,stats4]=regress(Y41',T4) % figure; subplot(3,2,2) plot(Ya,'k-') % Y42=Ym.*((ET1).^0.209.*(ET2).^0.7025.*(ET3).^0.2199.*(ET4).^0.1523) % hold on;plot(Y42,'b-') Y43=Ym.*((ET1).^b4(1).*(ET2).^b4(2).*(ET3).^b4(3).*(ET4).^b4(4)) hold on;plot(Y43,'r-') title('Jensen模型') xlabel('蒸发蒸腾量(mm)') ylabel('产量(kg/hm^2)') % legend('真实值','函数拟合值') %% 乘法模型2-Minhas模型 X51=log(1-(1-ET1).^2);X52=log(1-(1-ET2).^2);X53=log(1-(1-ET3).^2);X54=log(1-(1-ET4).^2); Y51=log(Ya/Ym); T5=[X51' X52' X53' X54']; [b5,bint5,r5,rint5,stats5]=regress(Y51',T5) % figure; subplot(3,2,4) plot(Ya,'k-') Y52=Ym.*((1-(1-ET1).^2).^b5(1).*(1-(1-ET2).^2).^b5(2).*(1-(1-ET3).^2).^b5(3).*(1-(1-ET1).^2).^b5(4)); hold on;plot(Y52,'r-') title('Minhas模型') xlabel('蒸发蒸腾量(mm)') ylabel('产量(kg/hm^2)') % legend('真实值','函数拟合值') %% 乘法模型3-Rao模型 X61=log(1-b2(1)*(1-ET1));X62=log(1-b2(2)*(1-ET2));X63=log(1-b2(3)*(1-ET3));X64=log(1-b2(4)*(1-ET4)); Y61=log(Ya/Ym); T6=[X61' X62' X63' X64']; [b6,bint6,r6,rint6,stats6]=regress(Y61',T6) % figure; subplot(3,2,6) plot(Ya,'k-') Y62=Ym.*((1-b2(1)*(1-ET1)).^b6(1).*(1-b2(2)*(1-ET2)).^b6(2).*(1-b2(3)*(1-ET3)).^b6(3).*(1-b2(4)*(1-ET4)).^b6(4)); hold on;plot(Y62,'r-') title('Rao模型') xlabel('蒸发蒸腾量(mm)') ylabel('产量(kg/hm^2)') % legend('真实值','函数拟合值') %% 平均绝对偏离和标准差 diff1=sum(abs(Y13-Ya))/11 diff2=sum(abs(Y22-Ya))/11 diff3=sum(abs(Y32-Ya))/11 diff4=sum(abs(Y43-Ya))/11 diff5=sum(abs(Y52-Ya))/11 diff6=sum(abs(Y62-Ya))/11 nor1=norm(abs(Y13-Ya)) nor2=norm(abs(Y22-Ya)) nor3=norm(abs(Y32-Ya)) nor4=norm(abs(Y43-Ya)) nor5=norm(abs(Y52-Ya)) nor6=norm(abs(Y62-Ya))