**********************************************************************;
* BIVARIATE.SAS -- Bivariate Normal Confidence Regions and Intervals *;
* Barry Moser, Dept. Experimental Statistics, LSU, Baton Rouge.      *;
* CREATED: 09/19/92                                                  *;
**********************************************************************;
Options PS=55 LS=78 PageNo=1 NoDate;
GOptions FText=SWISSX HText=1.0 CText=BLACK
         HTitle=1.0 FTitle=SWISSX CTitle=BLACK
         NoPrompt;

PROC IML WORKSIZE=100;
START CONTOUR(Contour,N,MEAN,COV,C,NPOINTS);
    /* find principle axes of ellipses */
    CALL EIGEN(V,E,(COV/N));
    /* set contour levels */
    A=SQRT(C*V[1]); B=SQRT(C*V[2]);
    /* parameterize the ellipse by angle */
    T=( (1:NPOINTS) - {1})#ATAN(1)#8/(NPOINTS-1);
    S=SIN(T); T=COS(T);
    S=S`*A; T=T`*B;
    /* form contour points */
    S =((E*(SHAPE(S,1)//SHAPE(T,1)))+MEAN@J(1,NPOINTS*NCOL(C),1) )`;
    Contour=SHAPE( S , NPOINTS ) ;
    /* returned as ncol pairs of columns for contours */
    FREE A B T S;
FINISH;

Start Interval(IntName,N,Mean,Var,C);
 Int=J(2,2,0);
 int[,1]=mean-sqrt(var/n)*C; /* lower limit */
 int[,2]=mean+sqrt(var/n)*C; /* upper limit */
 Print c, int;
 intc=RowCat(Char(Int,10,5));
 n1=rowcatc(IntName||{"X"});
 Call Symput(n1,intc[1]);
 n1=rowcatc(IntName||{"Y"});
 Call Symput(n1,intc[2]);
Finish;

Reset NoPrint Log;
N=42;
 /* CORRECT FOR THE MEAN */
MEAN={2.1905, 10.0476};
COV={1.1823 1.0883, 1.0883 11.3635};
Print Mean COV;
C=(2#(N-1)/(N-2))#FInv(0.95,2,(N-2)); /* F with 2,n-2 d.f. */
Run Contour(Results,N,MEAN,COV,C,100);
ColNames={"X" "Y"};
Create BC From Results[Colname=ColNames];
Append From Results;
Close BC;
Free ColNames;
 /* Now Compute Univariate Confidence Intervals */
var=vecdiag(Cov);
 /* T2 intervals */
C=sqrt((2#(N-1)/(N-2))#FInv(0.95,2,(N-2))); /* F with 2,n-2 d.f. */
Run Interval({"TSQ"},N,Mean,Var,C);
 /* Bonferroni Intervals */
C=TInv(1-0.05/4,(N-1));
Run Interval({"BON"},N,Mean,Var,C);
 /* Usual Univariate Intervals */
C=TInv(1-0.05/2,(N-1));
Run Interval({"UNI"},N,Mean,Var,C);
MeanC=Char(Mean,10,5);
Call Symput("MeanX",meanC[1]);
Call Symput("MeanY",meanC[2]);
Quit;

Data Ann;
 Drop SX SY;
 Length Text $10.;
 Size=1.5;
 Function="LABEL   ";
 HSYS="4"; YSYS="2"; XSYS="2";
 Position="B";
 Color="BLACK   ";
 Text="T-Square";
 SX="&TSQX"; SY="&TSQY";
 Link GetXY;
 Text="Bonferroni";
 SX="&BONX"; SY="&BONY";
 Link GetXY;
 Text="Univariate";
 SX="&UNIX"; SY="&UNIY";
 Link GetXY;
 X=input(substr("&MeanX",1,10),10.5);
 Y=input(substr("&MeanY",1,10),10.5);
 Text="m";
 Style="GREEK   ";
 Position="2";
 Output;
 Text="^";
 Style=' ';
 Position="2";
 Output;
 Function="SYMBOL";
 Text="DOT";
 Style=" ";
 Position="5";
 Size=1;
 Output;
 Stop;
GETXY:
 X=(input(substr(SX,1,10),10.5)+input(substr(SX,11,10),10.5))/2;
 Y=input(substr(SY,11,10),10.5);
 Output;
Return;
Run;

Title1 H=1.5 F=SWISSX 'Bivariate Normal 95% Confidence Ellipse';
Title2 H=1.5 F=SWISSX 'With Various 95% Univariate Confidence Intervals';
Proc GPlot Data=BC Annotate=ANN;
Plot Y*X=1 / HRef=&TSQX &BONX &UNIX VREF=&TSQY &BONY &UNIY;
Symbol1 c=Black i=join v=none r=1 h=1.5;
Run; Quit;

Title2 H=1.5 F=SWISSX 'With 95% T-Square Confidence Intervals';
Proc GPlot Data=BC;
Plot Y*X=1 / HRef=&TSQX VREF=&TSQY;
Symbol1 c=Black i=join v=none r=1 h=1.5;
Run; Quit;

Title2 H=1.5 F=SWISSX 'With 95% Bonferroni Confidence Intervals';
Proc GPlot Data=BC;
Plot Y*X=1 / HRef=&BONX VREF=&BONY;
Symbol1 c=Black i=join v=none r=1 h=1.5;
Run; Quit;

Title2 H=1.5 F=SWISSX 'With Usual 95% Univariate Confidence Intervals';
Proc GPlot Data=BC;
Plot Y*X=1 / HRef=&UNIX VREF=&UNIY;
Symbol1 c=Black i=join v=none r=1 h=1.5;
Run; Quit;