X=[ones(size(x1),1),x1,x2]; %stepwise(x,y);
[b1,bint1,r1,rint1,s1]=regress(y,X);%regress of cd4 after sort xx=-10:80;
yy=b1(1)+b1(2)*xx+b1(3)*xx.^2;yyn=max(yy);
xxn=(sqrt(b1(2)^2-4*b1(3)*(b1(1)-yyn))-b1(2))/(2*b1(3)); plot(xx,yy) max(yy)
%plot(data2(:,2),data2(:,3),'.');
hold on x3=data2(:,4); x4=data2(:,4).^2;
y2=data2(:,5)*50; %enlarge y2 to '50*y2' so it %will be clearer when ploted with yy; xs=[x3,x4];
Xs=[ones(size(x3),1),x3,x4]; %stepwise(x,y);
[b2,bint2,r2,rint2,s2]=regress(y2,Xs);%regress of hiv after sort xx=-10:80;
yy=b2(1)+b2(2)*xx+b2(3)*xx.^2;yyn=min(yy);
xxn2=(sqrt(b2(2)^2-4*b2(3)*(b2(1)-yyn))-b2(2))/(2*b2(3)); plot(xx,yy) min(yy) grid
text(60,320,'\\leftarrow HIV'); text(60,105,'\\leftarrow CD4'); title('HIV and CD4'); xlabel('date(week)'); ylabel('HIV*50 or CD4'); hold off
xt1=data2(:,3);xt2=xt1.^2;xt3=1./(xt1+1);yt=data2(:,5); xt=[xt1,xt2];
Xt=[ones(size(xt1),1),xt1,xt2,xt3]; %stepwise(xt,yt);
[b3,bint3,r3,rint3,s3]=regress(yt,Xt); xx=0.0000001:300;
yyt=b3(1)+b3(2)*xx+b3(3)*xx.^2+b3(4)./(xx+1);
figure(2);
plot(xx,yyt) hold on
21
plot(data2(:,3),data2(:,5),'.') hold off
xxn úte when cd4 reach max(cd4) xxn2 úte when hiv reach min(hiv)
figure(3);
bx1=data1(:,2);bx2=bx1.^2;by=data1(:,3); bX=[ones(size(bx1),1),bx1,bx2]; bx=[bx1,bx2]; %stepwise(bx,by)
[b4,bint4,r4,rint4,s4]=regress(by,bX);%regress of cd4 with no sort xx=-10:80;
yy1=b4(1)+b4(2)*xx+b4(3)*xx.^2;
bx1=data1(:,4);bx2=bx1.^2;by=data1(:,5)*50; bX=[ones(size(bx1),1),bx1,bx2]; bx=[bx1,bx2]; %stepwise(bx,by)
[b5,bint5,r5,rint5,s5]=regress(by,bX);%regress of hiv with no sort yy2=b5(1)+b5(2)*xx+b5(3)*xx.^2; %plot(xx,yy1,xx,yy2,bx1,by,'.') plot(xx,yy1,xx,yy2) grid
text(60,320,'\\leftarrow HIV'); text(61,103,'\\leftarrow CD4'); title('HIV and CD4'); xlabel('date(week)'); ylabel('HIV*50 or CD4'); yyn=max(yy1)
xxn3=(sqrt(b4(2)^2-4*b4(3)*(b4(1)-yyn))-b4(2))/(2*b4(3)); yyn=min(yy2)
xxn4=(sqrt(b5(2)^2-4*b5(3)*(b5(1)-yyn))-b5(2))/(2*b5(3)); xxn3 xxn4
8.1.2 问题二回归分析的程序代码
size(data201);
[ra,ir]=sort(data201(:,2)); D=[ra,ir]; for k=1:4 j(k)=0; for i=1:5036 if D(i,1)==k
22
j(k)=j(k)+1; end end end
X1=data201(D(1:j(1),2),:);
X2=data201(D(j(1)+1:j(1)+j(2),2),:);
X3=data201(D(j(1)+j(2)+1:j(1)+j(2)+j(3),2),:);
X4=data201(D(j(1)+j(2)+j(3)+1:j(1)+j(2)+j(3)+j(4),2),:);
data1=X4; n=size(data1);
[ra,ir]=sort(data1(:,2)); W=[ra,ir]; k=1; for i=1:n(1)
%if (data1(i,4)<=1&data1(i,5)>3&data1(i,5)<4.5 if (data1(i,4)<=1)&data1(i,5)<3 %if (data1(i,4)<=1&data1(i,5)>4.5 s(k)=i; k=k+1; end end s';
M=data1(s',:); j=0; for i=2:n(1)
if data1(i,1)-data1(i-1,1)==1 j=j+1; end end n2=size(M); k=1; for i=1:n2(1) for j=1:n(1)
if data1(j,1)==M(i,1) ss(k)=j; k=k+1; end end end ss;
23
data2=(data1(ss,:)); x1=data1(:,3); x2=x1.^2; x3=data1(:,4); x4=x3.^2; x5=x1.*x3; y=data1(:,5); xx=[x1,x2,x3,x4,x5]; x=[x1,x3];
rstool(x,y,'purequadratic');
8.1.3 问题二方差分析的程序代码
n=size(data201); k=1; for i=1:n(1)
if data201(i,4)>=22&data201(i,4)<=26 data20(k,:)=data201(i,:); k=k+1; end end data20; n2=size(data20);
%for j=1:4 k=1; for i=1:n2(1) if data20(i,2)==1
data21(k,:)=data20(i,:); k=k+1; end end X1=data21;
clear data21; k=1; for i=1:n2(1) if data20(i,2)==2
data21(k,:)=data20(i,:); k=k+1; end end X2=data21;
24