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