1 1 AN ADJUSTED LIKELIHOOD RATIO ALGORITHM FOR THE BEHRENS-FISHER PROBLEM HAMPARSUM BOZDOGAN AND DONALD E. RAMIREZ DEPARTMENT OF MATHEMATICS MATH-ASTRO BUILDING UNIVERSITY OF VIRGINIA CHARLOTTESVILLE, VIRGINIA 22903 U.S.A. TECHNICAL REPORT 18 (FEBRUARY 3,1987) 1 INTRODUCTION THIS TECHNICAL REPORT PRESENTS THE PROGRAM BFSOLV WHICH IS DEVELOPED IN OUR PAPER: AN ADJUSTED LIKELIHOOD ALGORITHM FOR THE BEHRENS-FISHER PROBLEM, STATISTISCHE HEFTE (TO APPEAR) BFSOLV IS A FORTRAN 77 PROGRAM WHICH FINDS THE MAXIMUM LIKELIHOOD ESTIMATOR FOR THE COMMON MEAN FOR RANDOM NORMAL DATA COMING FROM TWO GROUPS WITH COMMON MEANS BUT WITH DIFFERENT VARIANCES. THE OUTPUT IS THE SIGNIFICANCE LEVELS FOR THE WELCH PROCEDURE, FOR THE SCALED LIKELIHOOD PROCEDURE, AND FOR THE ADJUSTED LIKELIHOOD PROCEDURE. THE USER INPUTS TO THE SUBROUTINE BFSOLV ARE N1 THE NUMBER OF CASES FOR THE FIRST GROUP, N2 THE NUMBER OF CASES FOR THE SECOND GROUP, AND THE DATA X1 AND X2. INDEX APPENDIX A CONTAINS A DATA SET WITH N1 = 5 AND N2 = 20 APPENDIX B CONTAINS THE TERMINAL OUTPUT FROM BFSOLV APPENDIX C CONTAINS THE FULL OUTPUT FROM THE PROGRAM BFSOLV APPENDIX D CONTAINS THE FORTRAN PROGRAM BFSOLV APPENDIX E CONTAINS THE CALL STATEMENTS FROM BFSOLV A COPY OF THIS PROGRAM IS AVAILABLE ON A FLOPPY DISK FROM EITHER OF THE AUTHORS ON AN AS-IS-BASIS FOR $50.00. ************************************************************************ APPENDIX A: THE ORIGINAL DATA 5, 20 131.7 88.9 112.6 103.0 117.6 86.8 117.3 102.0 109.8 97.7 83.4 101.7 99.8 92.6 93.8 93.4 89.6 97.9 96.9 91.7 101.3 109.9 95.7 109.9 94.2 NOTE: N1 IS 5 AND X1 IS FROM N(MEAN=100,VARIANCE=300) N2 IS 20 AND X2 IS FROM N(MEAN=100,VARIANC= 60) ************************************************************************ APPENDIX B: TERMINAL OUTPUT FROM BFSOLV OK! SEG BF_FEB3 ENTER THE NAME OF THE OUTPUT FILE BF_OUT PROGRAM COMPUTES SIGNIFICANCE LEVEL FOR THE BEHRENS-FISHER PROBLEM ENTER NAME OF DATA FILE CONTAINING THE PAIR N1,N2 AND COLUMN OF DATA BF_DATA READING DATA DATA READ **** STOP ************************************************************************ APPENDIX C: FULL OUTPUT FROM BFSOLV OUTFILE = BF_OUT PROGRAM COMPUTES SIGNIFICANCE LEVEL FOR THE BEHRENS-FISHER PROBLEM DATA FILE = BF_DATA ************************************************************************ FILE HAS SAMPLE SIZES = 5 20 GMEAN1 GMEAN2 GVAR1 GVAR2 110.75998 98.269913 205.34937 68.764679 SAMPLE SIZES ARE 5 20 THE MLE ROOTS ARE UNIQUE FOR THE ALTERNATIVE HYPOTHESIS: MEAN1 MEAN2 VAR1 VAR2 110.76 98.270 205.35 68.765 FOR THE NULL HYPOTHESIS: MEAN VAR1 VAR2 98.869 346.74 69.117 FOR THE WELCH PROCEDURE: T-VALUE DF P-VALUE 1.6848261 4.57907 0.15820 FOR THE SCALED LIKELIHOOD RATIO TEST: -2*SCALED-LOGLIKELIHOOD DF P-VALUE 1.5339556 1.00000 0.21552 FOR THE LIKELIHOOD RATIO TEST ADJUSTED FOR THE FIRST MOMENT: -2*LOGLIKELIHOOD DF P-VALUE 2.7234950 1.77547 0.21754 PROGRAM TERMINATED NORMALLY ************************************************************************ APPENDIX D: THE FORTRAN PROGRAM BFSOLV PROGRAM MAIN C DEVELOPED BY HAMPARSUM BOZDOGAN AND DONALD E. RAMIREZ C MATHEMATICS DEPARTMENT, MATH-ASTRO BUILDING C UNIVERSITY OF VIRGINIA, CHARLOTTESVILLE, VA. 22903, U.S.A. C COPYRIGHT @ 02/03/87 INTEGER N1,N2,OUTUN,INUN,MAXNO PARAMETER (OUTUN=6,INUN=5,MAXNO=1000) REAL X1,X2,MEAN1,MEAN2,VAR1,VAR2, + BFMEAN(2),BFVAR1,BFVAR2,ALSCAL,BIAS, + TVALW,ALPHAW,DFW,M2LNLK,ALADJ DIMENSION X1(MAXNO),X2(MAXNO) LOGICAL UNIQUE,ILLCND EXTERNAL BFSOLV,READIN,PRTOUT,FLOPEN C FLOPEN OPENS A FILE WITH UNIT NUMBER OUTUN CALL FLOPEN(OUTUN) PRINT 700 WRITE(OUTUN,700) 700 FORMAT(/,1X,'PROGRAM COMPUTES SIGNIFICANCE LEVEL FOR THE ', + 'BEHRENS-FISHER PROBLEM') C READIN READS SAMPLE SIZES AND DATA FROM A FILE C FILE MUST CONTAIN SAMPLE SIZES ON THE FIRST LINE C FOLLOWED BY A SINGLE COLUMN OF DATA CALL READIN(INUN,OUTUN,N1,N2,X1,X2) C BFSOLV COMPUTES THE SIGNIFICANCE LEVELS FOR THE WELCH PROCEDURE C FOR THE SCALED LIKELIHOOD PROCEDURE AND C FOR THE ADJUSTED LIKELIHOOD PROCEDURE C N1 INPUT: THE SAMPLE SIZE FOR THE FIRST GROUP C N2 INPUT: THE SAMPLE SIZE FOR THE SECOND GROUP C X1 INPUT: THE VECTOR OF DATA FOR THE FIRST GROUP C X2 INPUT: THE VECTOR OF DATA FOR THE SECOND GROUP C MEAN1 OUTPUT: THE SAMPLE MEAN FOR THE FIRST GROUP C MEAN2 OUTPUT: THE SAMPLE MEAN FOR THE SECOND GROUP C VAR1 OUTPUT: THE BIASED VARIANCE FOR THE FIRST GROUP C VAR2 OUTPUT: THE BIASED VARIANCE FOR THE SECOND GROUP C BFMEAN OUTPUT: VECTOR OF LENGTH 2 FOR MLE FOR THE COMMON MEAN C BFVAR1 OUTPUT: THE MLE FOR THE VARIANCE FOR THE FIRST GROUP C BFVAR2 OUTPUT: THE MLE FOR THE VARIANCE FOR THE SECOND GROUP C M2LNLK OUTPUT: MINUS TWICE THE LOG-LIKELIHOOD RATIO C BIAS OUTPUT: EXPECTED VALUE OF M2LNLK C ALSCAL OUTPUT: THE SIGNIFICANCE LEVEL FOR THE SCALED LIKELIHOOD C PROCEDURE C TVALW OUTPUT: THE WELCH T-STATISTIC C DFW OUTPUT: THE DF'S FOR THE WELCH PROCEDURE C ALPHAW OUTPUT: THE SIGNIFICANCE LEVEL FOR THE WELCH PROCEDURE C ALADJ OUTPUT: THE SIGNIFICANCE LEVEL FOR THE ADJUSTED LIKELIHOOD C PROCEDURE C UNIQUE OUTPUT: LOGICAL VARIABLE FOR WHETHER B-F MEAN IS UNIQUE C ILLCND OUTPUT: LOGICAL VARIABLE FOR WHETHER THE B-F EQUATION C MAY BE ILL-CONDITIONED CALL BFSOLV(N1,N2,X1,X2,MEAN1,MEAN2,VAR1,VAR2,BFMEAN, + BFVAR1,BFVAR2,M2LNLK,BIAS,ALSCAL,TVALW,DFW,ALPHAW, + ALADJ,UNIQUE,ILLCND) C PRTOUT PRINTS STATISTICS FROM BFSOLVE CALL PRTOUT(N1,N2,MEAN1,MEAN2,VAR1,VAR2,BFMEAN,BFVAR1,BFVAR2, + M2LNLK,BIAS,ALSCAL,TVALW,DFW,ALPHAW,ALADJ,UNIQUE, + ILLCND,OUTUN) WRITE(OUTUN,710) 710 FORMAT(/1X,'PROGRAM TERMINATED NORMALLY') CLOSE(OUTUN) END ********************************************************************* SUBROUTINE FLOPEN(OUTUN) INTEGER OUTUN CHARACTER*40 OUTFIL LOGICAL EX 10 CONTINUE PRINT 700 700 FORMAT(/,1X,'ENTER THE NAME OF THE OUTPUT FILE') READ 710,OUTFIL 710 FORMAT(A) INQUIRE(FILE=OUTFIL,EXIST=EX) IF (EX) THEN PRINT 720,OUTFIL 720 FORMAT(/1X,'FILE = ',A,' EXISTS. TRY AGAIN.') GOTO 10 END IF OPEN(FILE=OUTFIL,UNIT=OUTUN,ERR=20) WRITE(OUTUN,730) OUTFIL 730 FORMAT(/,1X,'OUTFILE = ',A) RETURN 20 CONTINUE PRINT 740,OUTFIL 740 FORMAT(/,1X,'ERROR IN OPENING FILE = ',A) GOTO 10 END ********************************************************************* SUBROUTINE READIN(INUN,OUTUN,N1,N2,X1,X2) REAL X1(*),X2(*) INTEGER I,N1,N2,INUN,OUTUN CHARACTER*40 INFILE LOGICAL EX 10 CONTINUE PRINT 700 700 FORMAT(/1X,'ENTER NAME OF DATA FILE CONTAINING THE PAIR N1,N2 ', + 'AND COLUMN OF DATA') READ(*,'(A)')INFILE INQUIRE(FILE=INFILE,EXIST=EX) IF (.NOT.EX) THEN PRINT*,'FILE = ',INFILE PRINT*,'DOES NOT EXIST. PLEASE TRY AGAIN.' GOTO 10 ENDIF WRITE(OUTUN,710) INFILE 710 FORMAT(/1X,'DATA FILE = ',A) WRITE(OUTUN,720) 720 FORMAT(/1X,72('*')) OPEN(UNIT=INUN,FILE=INFILE,ERR=40) PRINT 730 730 FORMAT(/1X,'READING DATA') READ(INUN,*) N1,N2 DO 20,I=1,N1 READ(INUN,*) X1(I) 20 CONTINUE DO 30,I=1,N2 READ(INUN,*) X2(I) 30 CONTINUE PRINT 740 740 FORMAT(/1X,'DATA READ') WRITE(OUTUN,750) N1,N2 750 FORMAT(/1X,'FILE HAS SAMPLE SIZES = ',2I5) CLOSE(INUN) RETURN 40 CONTINUE PRINT 760,INFILE 760 FORMAT(/,1X,'ERROR IN OPENING FILE = ') GOTO 10 END ******************************************************************** SUBROUTINE PRTOUT(N1,N2,MEAN1,MEAN2,VAR1,VAR2,BFMEAN,BFVAR1, + BFVAR2,M2LNLK,BIAS,ALSCAL,TVALW,DFW,ALPHAW, + ALADJ,UNIQUE,ILLCND,OUTUN) INTEGER OUTUN,N1,N2 REAL MEAN1,MEAN2,VAR1,VAR2, + BFMEAN(2),BFVAR1,BFVAR2,M2LNLK,BIAS,ALSCAL, + TVALW,DFW,ALPHAW,ALADJ,ONE LOGICAL UNIQUE,ILLCND ONE=1.0 WRITE(OUTUN,700) 'GMEAN1','GMEAN2','GVAR1','GVAR2' 700 FORMAT(/5X,A,9X,A,10X,A,10X,A) WRITE(OUTUN,710) MEAN1,MEAN2,VAR1,VAR2 710 FORMAT(4G15.8) WRITE(OUTUN,720) N1,N2 720 FORMAT(/1X,'SAMPLE SIZES ARE',2I5) IF (UNIQUE) THEN WRITE(OUTUN,730) 730 FORMAT(/1X,'THE MLE ROOTS ARE UNIQUE') IF (ILLCND) THEN WRITE(OUTUN,735) 735 FORMAT(/1X,72('*'),/1X,'WARNING: THE MLE EQUATIONS ', + 'MAY BE ILL-CONDITIONED!',/1X,72('*')) END IF ELSE WRITE(OUTUN,740) BFMEAN(1),BFMEAN(2) 740 FORMAT(/1X,72('*'), + /1X,'WARNING: THE MLE ROOTS ARE NOT UNIQUE', + /1X,'THE TWO MEANS ARE ',G15.5,5X,G15.5, + /1X,'LARGER MEAN IS USED BELOW', + /1X,'ALL SIGNIFICANCE LEVELS SHOULD BE SMALL', + /1X,72('*')) ENDIF WRITE(OUTUN,750) 750 FORMAT(/1X,'FOR THE ALTERNATIVE HYPOTHESIS:') WRITE(OUTUN,760) 'MEAN1','MEAN2','VAR1','VAR2' 760 FORMAT(7X,A,11X,A,12X,A,12X,A) WRITE(OUTUN,770) MEAN1,MEAN2,VAR1,VAR2 770 FORMAT(4(1X,G15.5)) WRITE(OUTUN,775) 775 FORMAT(/1X,'FOR THE NULL HYPOTHESIS:') WRITE(OUTUN,780) 'MEAN','VAR1','VAR2' 780 FORMAT(15X,A,21X,A,12X,A) WRITE(OUTUN,790) BFMEAN(1),BFVAR1,BFVAR2 790 FORMAT(8X,G15.5,9X,2(1X,G15.5)) WRITE(OUTUN,800) 800 FORMAT(/1X,'FOR THE WELCH PROCEDURE:') WRITE(OUTUN,810) 'T-VALUE','DF','P-VALUE' 810 FORMAT(7X,A,17X,A,18X,A) WRITE(OUTUN,820) TVALW,DFW,ALPHAW 820 FORMAT(G18.8,5X,F10.5,17X,F8.5) WRITE(OUTUN,830) 830 FORMAT(/1X,'FOR THE SCALED LIKELIHOOD RATIO TEST:') WRITE(OUTUN,840) '-2*SCALED-LOGLIKELIHOOD','DF','P-VALUE' 840 FORMAT(1X,A,7X,A,18X,A) WRITE(OUTUN,850) M2LNLK/BIAS, ONE, ALSCAL 850 FORMAT(G18.8,5X,F10.5,17X,F8.5) WRITE(OUTUN,855) 855 FORMAT(/1X,'FOR THE LIKELIHOOD RATIO TEST ADJUSTED FOR THE', + ' FIRST MOMENT:') WRITE(OUTUN,860) '-2*LOGLIKELIHOOD','DF','P-VALUE' 860 FORMAT(1X,A,14X,A,18X,A) WRITE(OUTUN,850) M2LNLK,BIAS,ALADJ C850 FORMAT(G18.8,5X,F10.5,17X,F8.5) PRINT 870 870 FORMAT(/1X) END ********************************************************************* SUBROUTINE BFSOLV(N1,N2,X1,X2,MEAN1,MEAN2,VAR1,VAR2,BFMEAN, + BFVAR1,BFVAR2,M2LNLK,BIAS,ALSCAL,TVALW,DFW,ALPHAW, + ALADJ,UNIQUE,ILLCND) C DEVELOPED BY HAMPARSUM BOZDOGAN AND DONALD E. RAMIREZ C MATHEMATICS DEPARTMENT, MATH-ASTRO BUILDING C UNIVERSITY OF VIRGINIA, CHARLOTTESVILLE, VA. 22903, U.S.A. C COPYRIGHT @ 02/03/87 C BFSOLV COMPUTES THE SIGNIFICANCE LEVELS FOR THE WELCH PROCEDURE C FOR THE SCALED LIKELIHOOD PROCEDURE AND C FOR THE ADJUSTED LIKELIHOOD PROCEDURE C N1 INPUT: THE SAMPLE SIZE FOR THE FIRST GROUP C N2 INPUT: THE SAMPLE SIZE FOR THE SECOND GROUP C X1 INPUT: THE VECTOR OF DATA FOR THE FIRST GROUP C X2 INPUT: THE VECTOR OF DATA FOR THE SECOND GROUP C MEAN1 OUTPUT: THE SAMPLE MEAN FOR THE FIRST GROUP C MEAN2 OUTPUT: THE SAMPLE MEAN FOR THE SECOND GROUP C VAR1 OUTPUT: THE BIASED VARIANCE FOR THE FIRST GROUP C VAR2 OUTPUT: THE BIASED VARIANCE FOR THE SECOND GROUP C BFMEAN OUTPUT: VECTOR OF LENGTH 2 FOR THE MLE FOR THE COMMON MEAN C BFVAR1 OUTPUT: THE MLE FOR THE VARIANCE FOR THE FIRST GROUP C BFVAR2 OUTPUT: THE MLE FOR THE VARIANCE FOR THE SECOND GROUP C M2LNLK OUTPUT: MINUS TWICE THE LOG-LIKELIHOOD RATIO C BIAS OUTPUT: EXPECTED VALUE OF M2LNLK C ALSCAL OUTPUT: THE SIGNIFICANCE LEVEL FOR THE SCALED LIKELIHOOD C PROCEDURE C TVALW OUTPUT: THE WELCH T-STATISTIC C DFW OUTPUT: THE DF'S FOR THE WELCH PROCEDURE C ALPHAW OUTPUT: THE SIGNIFICANCE LEVEL FOR THE WELCH PROCEDURE C ALADJ OUTPUT: THE SIGNIFICANCE LEVEL FOR THE ADJUSTED LIKELIHOOD C PROCEDURE C UNIQUE OUTPUT: LOGICAL VARIABLE FOR WHETHER BFMEAN IS UNIQUE C ILLCND OUTPUT: LOGICAL VARIABLE FOR WHETHER THE B-F EQUATION C MAY BE ILL-CONDITIONED INTEGER N1,N2,IER, + LEVOLD,LEVNEW REAL X1,X2,MEAN1,MEAN2,VAR1,VAR2, + BFMEAN(2),BFVAR1,BFVAR2,ALSCAL, + TWELCH,TVALW,ALPHAW,DFW, + M2LNLK,ALADJ,BIAS,ONE DIMENSION X1(*),X2(*) C UERSET,MDTD,MDCH ARE IMSL ROUTINES C UERSET: SETS MESSAGE LEVEL FOR IMSL ROUTINES C MDTD: COMPUTES TWO-SIDED TAIL REGIONS FOR THE T-DISTRIBUTION C MDCH: COMPUTES LEFT-HANDED TAIL REGIONS FOR THE CHI-SQUARE C DISTRIBUTION EXTERNAL STATS,UERSET,BFVAR,MDTD, + MDCH,MLUNIQ,TWELCH LOGICAL UNIQUE,ILLCND INTRINSIC LOG C STATS: COMPUTES MEAN AND BIASED VARIANCE CALL STATS(X1,N1,MEAN1,VAR1) CALL STATS(X2,N2,MEAN2,VAR2) LEVNEW = 1 CALL UERSET(LEVNEW,LEVOLD) C MLUNIQ: COMPUTES MLE FOR THE COMMON MEAN CALL MLUNIQ(UNIQUE,N1,N2,MEAN1,MEAN2,VAR1,VAR2,BFMEAN,ILLCND) C BFVAR: COMPUTES MLE'S FOR SEPARATE VARIANCES CALL BFVAR(X1,N1,X2,N2,BFMEAN(1),BFVAR1,BFVAR2) M2LNLK=N1*LOG(1.0+(BFMEAN(1)-MEAN1)**2/VAR1) + +N2*LOG(1.0+(BFMEAN(1)-MEAN2)**2/VAR2) TVALW=TWELCH(MEAN1,MEAN2,VAR1,VAR2,N1,N2,DFW) C MDTD: COMPUTES TWO-SIDED TAIL REGIONS FOR THE T-DISTRIBUTION CALL MDTD(TVALW,DFW,ALPHAW,IER) BIAS=DFW/(DFW-2.0) C MDCH: COMPUTES LEFT-HANDED TAIL REGIONS FOR THE CHI-SQUARE C DISTRIBUTION CALL MDCH(M2LNLK,BIAS,ALADJ,IER) ALADJ=1.0-ALADJ ONE=1.0 C MDCH: COMPUTES LEFT-HANDED TAIL REGIONS FOR THE CHI-SQUARE C DISTRIBUTION CALL MDCH(M2LNLK/BIAS,ONE,ALSCAL,IER) ALSCAL=1.0-ALSCAL RETURN END ****************************************************************** SUBROUTINE STATS(A,N,MEAN,VAR) REAL A,MEAN,VAR,SUM,SUMSQ INTEGER N,ROW DIMENSION A(*) INTRINSIC REAL SUM=0.0 SUMSQ=0.0 DO 10, ROW=1,N SUM=SUM+A(ROW) SUMSQ=SUMSQ+A(ROW)*A(ROW) 10 CONTINUE MEAN=SUM/REAL(N) VAR=(SUMSQ-SUM*SUM/REAL(N))/REAL(N) RETURN END ************************************************************************ SUBROUTINE MLUNIQ(UNIQUE,N1,N2,MEAN1,MEAN2,VAR1,VAR2,BFMEAN, + ILLCND) INTEGER N1,N2 REAL MEAN1,MEAN2,VAR1,VAR2,A,B,C,D,DISC REAL P,Q,R,AA,BB,SMALL,BFMEAN(2),W,Z,PHI,PI REAL LNLIKE,ROOT1,ROOT2,ROOT3,LTROOT,RTROOT,LTLNLK,RTLNLK REAL AADIV3,BBDIV2 INTRINSIC ABS,SQRT,SIGN,ACOS,COS,MIN,MAX EXTERNAL LNLIKE LOGICAL UNIQUE,ILLCND DATA SMALL /1E-7/ PI=ACOS(-1.0) C CUBIC IS A*X**3 + B*X**2 + C*X + D A=1.0/N1 + 1.0/N2 B=-((2.0/N1 + 1.0/N2)*MEAN1 + (1.0/N1 + 2.0/N2)*MEAN2) C=VAR1/N1 + VAR2/N2 + MEAN1*MEAN1/N1 + MEAN2*MEAN2/N2 + + 2.0*(1.0/N1 + 1.0/N2)*MEAN1*MEAN2 D=-(MEAN1*MEAN1*MEAN2/N1 + MEAN1*MEAN2*MEAN2/N2 + + VAR1*MEAN2/N1 + VAR2*MEAN1/N2) C TESTING FOR DISCRIMINANT=0 C CUBIC REWRITTEN AS X**3 + P*X**2 + Q*X + R P=B/A Q=C/A R=D/A AA=(3.0*Q-P*P)/3.0 AADIV3=AA/3.0 BB=(2.0*P*P*P-9.0*P*Q+27.0*R)/27.0 BBDIV2=BB/2.0 DISC=BB*BB/4.0 +AA*AA*AA/27.0 C REAL ROOTS OF A CUBIC ARE UNIQUE IF AND ONLY IF DISC > 0 C OR IF DISC = AA = BB = 0 C REFERENCE: HODGMAN, CHARLES D. (EDITOR), STANDARD MATHEMATICAL C TABLES (ELEVENTH EDITION), CHEMICAL RUBBER PUBLISHING C CO., 344-345 (1957). IF (DISC .GT. 0.0) THEN UNIQUE=.TRUE. ILLCND=.FALSE. W=(-BBDIV2+SQRT(DISC)) Z=(-BBDIV2-SQRT(DISC)) BFMEAN(1)=SIGN(1.0,W)*(ABS(W))**(1.0/3.0) + + SIGN(1.0,Z)*(ABS(Z))**(1.0/3.0) - P/3.0 BFMEAN(2)=BFMEAN(1) ELSE ILLCND=.TRUE. IF ((ABS(DISC) .LT. SMALL) .AND. (ABS(AA) .LT. SMALL) .AND. + (ABS(BB) .LT. SMALL)) THEN UNIQUE=.TRUE. BFMEAN(1)=-P/3.0 BFMEAN(2)=BFMEAN(1) ELSE PHI=ACOS(-BBDIV2/SQRT(-AA**3/27.0)) ROOT1=2.0*SQRT(-AADIV3)*COS(PHI/3.0) ROOT2=2.0*SQRT(-AADIV3)*COS((PHI+2.0*PI)/3.0) ROOT3=2.0*SQRT(-AADIV3)*COS((PHI+4.0*PI)/3.0) C POSSIBLE SOLUTIONS FOR THE MLE FOR THE COMMON MEAN C ARE MIN AND MAX OF THE THREE REAL ROOTS LTROOT=MIN(ROOT1,ROOT2,ROOT3) - P/3.0 RTROOT=MAX(ROOT1,ROOT2,ROOT3) - P/3.0 LTLNLK=LNLIKE(LTROOT,N1,N2,MEAN1,MEAN2,VAR1,VAR2) RTLNLK=LNLIKE(RTROOT,N1,N2,MEAN1,MEAN2,VAR1,VAR2) C IF LTLNLK=RTLNLK THEN BFMEAN(1) WILL BE SET TO RTROOT C AND BFMEAN(2) WILL BE LTROOT IF (LTLNLK.GT.RTLNLK) THEN BFMEAN(1)=LTROOT BFMEAN(2)=RTROOT ELSE BFMEAN(1)=RTROOT BFMEAN(2)=LTROOT ENDIF IF (ABS(LTLNLK-RTLNLK).GT.SMALL) THEN UNIQUE = .TRUE. ELSE UNIQUE = .FALSE. ENDIF ENDIF END IF RETURN END ************************************************************************ FUNCTION LNLIKE(MU,N1,N2,MEAN1,MEAN2,VAR1,VAR2) REAL LNLIKE,MU,MEAN1,MEAN2,VAR1,VAR2 INTEGER N1,N2 INTRINSIC LOG LNLIKE = -N1*LOG(VAR1 + (MEAN1-MU)**2) -N2*LOG(VAR2+(MEAN2-MU)**2) RETURN END ************************************************************************ SUBROUTINE BFVAR(Y1,K1,Y2,K2,BFMEAN,VAR1,VAR2) REAL Y1,Y2,BFMEAN,VAR1,VAR2,SUM INTEGER K1,K2,I DIMENSION Y1(*),Y2(*) INTRINSIC REAL SUM=0.0 DO 10, I=1,K1 SUM=SUM+(Y1(I)-BFMEAN)**2 10 CONTINUE VAR1=SUM/REAL(K1) SUM=0.0 DO 20, I=1,K2 SUM=SUM+(Y2(I)-BFMEAN)**2 20 CONTINUE VAR2=SUM/REAL(K2) RETURN END *********************************************************************** FUNCTION TWELCH(MEAN1,MEAN2,VAR1,VAR2,K1,K2,DFW) REAL TWELCH,MEAN1,MEAN2,VAR1,VAR2,UVAR1,UVAR2,DFW,A,B INTEGER K1,K2 INTRINSIC SQRT C UVAR1: THE UNBIASED VARIANCE FOR THE FIRST GROUP C UVAR2: THE UNBIASED VARIANCE FOR THE SECOND GROUP UVAR1=K1*VAR1/(K1-1) UVAR2=K2*VAR2/(K2-1) A=UVAR1/K1 B=UVAR2/K2 TWELCH=(MEAN1-MEAN2)/SQRT(A+B) DFW=(A+B)**2 DFW=DFW/(A*A/(K1-1)+B*B/(K2-1)) RETURN ********************************************************************** C SUBROUTINE MDTD(TVAL,Q,IER) C IMSL ROUTINE FOR TWO-SIDED STUDENT'S T PROBABILITIES C END ********************************************************************** C SUBROUTINE MDCH(CS,DF,P,IER) C IMSL ROUTINE FOR LEFT-HANDED CHI-SQUARED PROBABILITIES C END *********************************************************************** C SUBROUTINE UERSET(LEVNEW,LEVOLD) C SETS MESSAGE LEVEL FOR IMSL ROUTINES C END *********************************************************************** END *********************************************************************** APPENDIX E: THE CALL STATEMENTS FROM BFSOLVE 1: PROGRAM MAIN 15: CALL FLOPEN(OUTUN) 23: CALL READIN(INUN,OUTUN,N1,N2,X1,X2) 50: CALL BFSOLV(N1,N2,X1,X2,MEAN1,MEAN2,VAR1,VAR2,BFMEAN, 51: + BFVAR1,BFVAR2,M2LNLK,BIAS,ALSCAL,TVALW,DFW,ALPHAW, 52: + ALADJ,UNIQUE,ILLCND) 54: CALL PRTOUT(N1,N2,MEAN1,MEAN2,VAR1,VAR2,BFMEAN,BFVAR1,BFVAR2, 55: + M2LNLK,BIAS,ALSCAL,TVALW,DFW,ALPHAW,ALADJ,UNIQUE, 56: + ILLCND,OUTUN) 62: SUBROUTINE FLOPEN(OUTUN) 87: SUBROUTINE READIN(INUN,OUTUN,N1,N2,X1,X2) 129: SUBROUTINE PRTOUT(N1,N2,MEAN1,MEAN2,VAR1,VAR2,BFMEAN,BFVAR1, 130: + BFVAR2,M2LNLK,BIAS,ALSCAL,TVALW,DFW,ALPHAW, 131: + ALADJ,UNIQUE,ILLCND,OUTUN) 196: SUBROUTINE BFSOLV(N1,N2,X1,X2,MEAN1,MEAN2,VAR1,VAR2,BFMEAN, 197: + BFVAR1,BFVAR2,M2LNLK,BIAS,ALSCAL,TVALW,DFW,ALPHAW, 198: + ALADJ,UNIQUE,ILLCND) 246: CALL STATS(X1,N1,MEAN1,VAR1) 247: CALL STATS(X2,N2,MEAN2,VAR2) 249: CALL UERSET(LEVNEW,LEVOLD) 251: CALL MLUNIQ(UNIQUE,N1,N2,MEAN1,MEAN2,VAR1,VAR2,BFMEAN,ILLCND) 253: CALL BFVAR(X1,N1,X2,N2,BFMEAN(1),BFVAR1,BFVAR2) 258: CALL MDTD(TVALW,DFW,ALPHAW,IER) 262: CALL MDCH(M2LNLK,BIAS,ALADJ,IER) 267: CALL MDCH(M2LNLK/BIAS,ONE,ALSCAL,IER) 272: SUBROUTINE STATS(A,N,MEAN,VAR) 288: SUBROUTINE MLUNIQ(UNIQUE,N1,N2,MEAN1,MEAN2,VAR1,VAR2,BFMEAN, 289: + ILLCND) 375: SUBROUTINE BFVAR(Y1,K1,Y2,K2,BFMEAN,VAR1,VAR2) The following subroutines are from the IMSL package: UERSET MDTD MDCH ************************************************************************