function[beta,stats,ynew,ylr]=regres2(X,y,Xnew,pp) //定义函数 regres2，里面有四个变量:X(实 验各个因素的数据也可以称作观测值)；y(如实验测得的数据 F)；Xnew(新观测值,要新观测 的各个因素的值，如进给速度等)；pp(进行置信概率为 pp 的区间预测，一般取 0.95) beta=[];stats=[];ynew=[];ylr=[]; [n,p]=size(X); m=p+1; if n

=1) pp=0; end 者 pp 小于 0 或者 pp>=1，pp 置 0； X=log10(X); y=log10(y); Xnew=log10(Xnew); A=[ones(size(y)),X]; //如果参数的数目小于 4 或者 pp 不是数字或 //进行方差分析 //对输入的观测值和新观测值取 lg //ones(size(y))得到一个行数列数和 y 一样并且元 

//regress 用最小二乘法进行多元线性回归得到 beta //mean(y):y 的均值，SSR:回归平方和，自由度为 //定义 alpha (α) //定义模型 yhat=A*beta 素全部为 1 的矩阵，所以 A 为 X 和一列 1 的组成的矩阵 [beta,btml,rtm,rtml,stat]=regress(y,A); alpha=[0.05,0.01]; yhat=A*beta; SSR=(yhat-mean(y))'*(yhat-mean(y)); m-1 SSE=(yhat-y)'*(yhat-y); SST=(y-mean(y))'*(y-mean(y)); Fb=SSR/(m-1)/SSE*(n-m); Falpha=finv(1-alpha,m-1,n-m); 查表得到相应 Fa 的值，用于和 Fb 比较大小 table=cell(p+4,7); table(1,:)={'方差来源','偏差平方和','自由度','方差','F 值','Fa','显著性'};//各个列输出的内容 table(2+p,1:6)={'回归',SSR,m-1,SSR/(m-1),Fb,min(Falpha)}; table(3+p,1:6)={'剩余',SSE,n-m,SSE/(n-m),[],max(Falpha)}; table(4+p,1:3)={'总和',SST,n-1}; if Fb>max(Falpha) 显著’ //SSE:残差平方和，自由度为 n-m //SST:总离差平方和 ，自由度为 n-1 //定义 p+4 行，7 列的表格 //求得 F 值 //finv 为 F 分布的累加分布函数逆函数，相当于 //若 Fb 值大于 max(Falpha)，第七列输出‘高度 //若 F 值在 max(Falpha)和 min(Falpha)之间输出 //若小于 min(Falpha)或者大于 max(Falpha) table{2+p,7}='高度显著'; elseif (Fb<=max(Falpha))&(Fb>min(Falpha)) ‘显著’ table{2+p,7}='显著'; else 输出‘不显著’ table{2+p,7}='不显著'; end R2=SSR/SST; R=sqrt(R2); Sy=sqrt(SSE/(n-m)); mnX=mean(X); MNX=repmat(mnX,n,1); Ljj=diag((X-MNX)'*(X-MNX)); Pj=abs(beta(2:end).*sqrt(Ljj/SST)); c=diag(inv(A'*A));bj2=beta.*beta; SSj=bj2(2:end)./c(2:end); Fj=SSj/SSE*(n-m); Falpha=finv(1-[0.05,0.01],1,n-m); //判断结束 //进行相关系数检验 //对 R2 开平方得到 R(相关系数) //求得 Sy(剩余标准差),判断回归方程的精度 //进行回归系数的检验,即检验每个 Xi 对 y 的影响程度 //c 为矩阵 Xt*X 上的主对角元素 //偏回归平方和，越大说明 Xi 在回归方程中越重要 //求得 Fj 用于和 Falpha 的两个极限值 1,2 比较 //相当于查分布表得到 Falpha 的两个极限值 1,2 ind2=find(Fj>=Falpha(2)); ind1=find((Fj>=Falpha(1))&(Fj

xxx(ind1)=1; [tmp,zbx]=min(Fj); xzh={'不显著','显著','高度显著'}; for kk=1:p //定义 p 行的 table 要显示的内容 table(kk+1,:)={['x',num2str(kk)],SSj(kk),1,SSj(kk),Fj(kk),[],xzh{1+xxx(kk)}}; end table{2,6}=Falpha(1);table{3,6}=Falpha(2); stats={table,R,Sy,Pj}; //输出 table，R，Sy，Pj 的数值，判定每个 Xi 对 y 的影响程度 结束 if(nargin>2)&(isnumeric(Xnew)) //若输入参数数目大于二，并且 Xnew 是数字 [n1,p1]=size(Xnew); ynew 的置信区间 Xnew=[ones(n1,1),Xnew]; ynew=Xnew*beta; Shat2=SSE/(n-m)*(1+Xnew*inv(A'*A)*Xnew'); Syhat=sqrt(diag(Shat2)); //通过要观测的 Xnew 预计出 ta=tinv(0.5+pp/2,n-p-1); //查表得到在相应 pp 和自由度下的 t 分布的值 yl=ynew-ta*Syhat; yr=ynew+ta.*Syhat; ylr=[yl(:),yr(:)]; end //得到 ynew 的置信区间(yl--yr)