PROGRAM LL2CMP C C ESTIMATION OF TWO-COMPONENT ERROR MODEL FOR LINEAR CALIBRATION C FUNCTION OF THE FORM C C X = ALPHA + BETP MU EXP(ETA) + EPS C C WHERE MU IS THE TRUE CONCENTRATION AND ETA AND EPS ARE ERRORS C C THIS IS THE LATEST VERSION AS OF 9/10/97--BPD C THIS IS THE LATEST VERSION AS OF 12/2/97--DMR C THIS IS THE LATEST VERSION AS OF 7/9/98--MW AND BPD C THIS IS THE LATEST VERSION AS OF 8/19/98--BPD C C VERSION WITHOUT DIALOG BOX--BPD 2/23/99 USE MSIMSL IMPLICIT NONE REAL*8 ALPHA,BETP,SETA,SEPS,X(1000),MU(1000) REAL*8 LLM, C INTEGER N,NPTS,I, IPERM(1000) CHARACTER*12 INFILE,OUTFILE COMMON /DATA/ X,MU,N,NPTS PRINT *, 'ENTER INPUT FILE NAME' READ *, INFILE PRINT *, 'ENTER OUTPUT FILE NAME' READ *,OUTFILE OPEN (20,FILE=INFILE,STATUS='OLD') OPEN (23,FILE=OUTFILE) DO 7 I=1,1000 READ (20,*,END=5) MU(I),X(I) 7 CONTINUE 5 CONTINUE NPTS=18 N=I-1 DO 8 I=1,N IPERM(I)=I 8 CONTINUE CALL DSVRGP(N,MU,MU,IPERM) CALL DPERMU(N,X,IPERM,1,X) C=MU(N)/X(N) DO 23 I = 1, N X(I) = C*X(I) 23 CONTINUE CALL STARTVALS(I,MU, X, N,ALPHA,BETP,SETA,SEPS) IF (N.LE.0) GOTO 20 CALL LLMAX(ALPHA,BETP,SETA,SEPS,LLM) ALPHA = ALPHA/C BETP = BETP/C SEPS = SEPS/C WRITE(23,1000)ALPHA,BETP,SETA,SEPS,LLM 1000 FORMAT(' ALPHA = ',G17.10/ 1 ' BETP = ',G17.10/ 1 ' SIGMA_ETA = ',G17.10/ 1 ' SIGMA_EPSILON = ',G17.10/ 1 ' MAXIMIZED LOG LIKELIHOOD = ',G17.10) 20 STOP END SUBROUTINE LLMAX(ALPHA,BETP,SETA,SEPS,LLM) IMPLICIT NONE REAL*8 ALPHA,BETP,SETA,SEPS,LLM EXTERNAL LLTOT REAL*8 THETA0(4),S,FTOL,X(4),FVALUE INTEGER MAXFCN THETA0(1)=ALPHA THETA0(2)=BETP THETA0(3)=LOG(SETA) THETA0(4)=LOG(SEPS) S=.1 FTOL=1.0E-5 MAXFCN=1000 C DUMPOL IS AN IMSL NUMERICAL OPTIMIZATION ROUTINE THAT C USES POLYTOPE METHODS CALL DUMPOL(LLTOT,4,THETA0,S,FTOL,MAXFCN,X,FVALUE) ALPHA=X(1) BETP=X(2) SETA=EXP(X(3)) SEPS=EXP(X(4)) LLM=FVALUE RETURN END SUBROUTINE LLTOT(NPAR,THETA,LLT) IMPLICIT NONE REAL*8 THETA(4),LLT REAL*8 ALPHA,BETP,SETA,SEPS,TEMP,LL,X0,MU0 REAL*8 X(1000),MU(1000) INTEGER N,NPTS,NPTS0,NPAR,I COMMON /DATA/ X,MU,N,NPTS ALPHA=THETA(1) BETP=THETA(2) SETA=EXP(THETA(3)) SEPS=EXP(THETA(4)) NPTS0=NPTS TEMP=0 DO 10 I=1,N X0=X(I) MU0=MU(I) TEMP=TEMP+LL(X0,ALPHA,BETP,MU0,SETA,SEPS,NPTS0) 10 CONTINUE LLT=-TEMP RETURN END REAL*8 FUNCTION LL(X,ALPHA,BETP,MU,SETA,SEPS,NPTS) IMPLICIT NONE REAL*8 X,ALPHA,BETP,MU,SETA,SEPS REAL*8 A,B,QXFIX(1),QX(1000),QW(1000) REAL*8 PI,LLTEMP,X0,W0 REAL*8 F,FP,FPP,ETA,U INTEGER NPTS,IWEIGH,NFIX,J PI=3.1415926535 CALL NSOLV(X,ALPHA,BETP,MU,SETA,SEPS,ETA,F,FP,FPP) NFIX=0 IWEIGH=4 A=0 B=0 C DGQRUL IS AN IMSL ROUTINE FOR EVALUATING INTEGRALS C USING GAUSS-HERMITE QUADRATURE METHODS CALL DGQRUL(NPTS,IWEIGH,A,B,NFIX,QXFIX,QX,QW) LL=0.0 DO 10 J=1,NPTS U=QX(J) X0=U*SQRT(2./FPP)+ETA W0=QW(J) LLTEMP=-(X-ALPHA-BETP*MU*EXP(X0))**2/(2*SEPS**2) LLTEMP=LLTEMP-0.5*X0**2/SETA**2 LLTEMP=EXP(LLTEMP) LLTEMP=1/(2*PI*SEPS*SETA)*SQRT(2./FPP)*LLTEMP*W0*exp(U**2) LL=LL+LLTEMP 10 CONTINUE IF(LL.LE.0.0D0) THEN LL=-1.0D10 ELSE LL=LOG(LL) ENDIF RETURN END SUBROUTINE NSOLV(X,ALPHA,BETP,MU,SETA,SEPS,ETA,F,FP,FPP) IMPLICIT NONE REAL*8 X,ALPHA,BETP,MU,SETA,SEPS,PI REAL*8 ETA,F,FP,FPP,A,B,ERRREL,GTOL REAL*8 ETAOUT,FUN1,FUN1P REAL*8 X1,ALPHA1,BETP1,MU1,SETA1,SEPS1,FUN1OUT,FUN1POUT INTEGER MAXFN EXTERNAL FUN1,FUN1P COMMON /CONSTANTS/ X1,ALPHA1,BETP1,MU1,SETA1,SEPS1 C X1=X ALPHA1=ALPHA BETP1=BETP MU1=MU SETA1=SETA SEPS1=SEPS ERRREL=.001 PI = 3.1415926535 A=-30*seta B=30*seta ETA=(A+B)/2 GTOL=10D-8 MAXFN=500 !MAX NUM FUNEVALS C DUVMID IS AN IMSL ROUTINE FOR NUMERICAL EVALUATION OF C DERIVATIVES CALL DUVMID(FUN1,FUN1P,ETA,ERRREL,GTOL,MAXFN,A,B,ETAOUT,FUN1OUT, ! FUN1POUT) ETA=ETAOUT F=-log(2*PI*SEPS*SETA)+0.5*ETA**2/SETA**2+0.5*(X-ALPHA-BETP*MU* ! EXP(ETA))**2/SEPS**2 FP=-ETA/SETA**2-BETP*MU*EXP(ETA)*(X-ALPHA-BETP*MU*EXP(ETA))/ 1 SEPS**2 FPP=1/SETA**2+BETP**2*MU**2*EXP(2*ETA)/SEPS**2 FPP=FPP-BETP*MU*EXP(ETA)*(X-ALPHA-BETP*MU*EXP(ETA))/SEPS**2 RETURN END C EVALUATING FUCTION FOR MAXIMIZATION OF INTEGRAND REAL*8 FUNCTION FUN1(ETA) REAL*8 ETA,X1,ALPHA1,BETP1,MU1,SETA1,SEPS1,pi COMMON /CONSTANTS/ X1,ALPHA1,BETP1,MU1,SETA1,SEPS1 PI = 3.1415926535 FUN1=-log(2*pi*seps1*seta1)+0.5*ETA**2/SETA1**2+0.5*(X1-ALPHA1- ! BETP1*MU1*EXP(ETA))**2/SEPS1**2 RETURN END C EVALUATING DERIVATIVE FOR MAXIMIZATION OF INTEGRAND REAL*8 FUNCTION FUN1P(ETA) REAL*8 ETA,X1,ALPHA1,BETP1,MU1,SETA1,SEPS1 COMMON /CONSTANTS/ X1,ALPHA1,BETP1,MU1,SETA1,SEPS1 FUN1P=ETA/SETA1**2-BETP1*MU1*EXP(ETA)*(X1-ALPHA1-BETP1*MU1*EXP(ETA 1 ))/SEPS1**2 RETURN END C STARTVALS GENERATES APPROPRIATE STARTING PARAMETERS C FROM THE DATA USING WEIGHTED LEAST SQUARES SUBROUTINE STARTVALS(I,MU,X,N,ALPHA,BETP,SETA,SEPS) USE MSIMSL IMPLICIT NONE REAL*8 MU(1000),X(1000),Y(1000),Z(1000),LOGARRAY(1000),SEPS REAL*8 LOWMU,HIGHMU,SETA,ALPHA,BETP INTEGER I,J,K,L,N,COUNT LOWMU=MU(1) Y(1)=X(1) J=1 N=I-1 5 CONTINUE 10 IF (MU(J+1).EQ.LOWMU) THEN Y(J+1)=X(J+1) J=J+1 GOTO 10 ENDIF IF (J.LE.1) THEN SEPS=1.0 GOTO 30 ENDIF CALL STANDDEV(Y,J,SEPS) 30 CONTINUE C FINDING ESTIMATE FOR SETA N=I-1 HIGHMU=MU(N) Z(1)=X(N) K=N COUNT=1 20 IF (MU(K-1).EQ.HIGHMU) THEN Z(COUNT+1)=X(K-1) K=K-1 COUNT=COUNT+1 GOTO 20 ENDIF IF (COUNT.LE.1) THEN SETA=0.4 GOTO 35 ENDIF DO 25 L=1,COUNT LOGARRAY(L)=LOG(Z(L)) 25 CONTINUE CALL STANDDEV(LOGARRAY,COUNT,SETA) 35 CONTINUE CALL WLEASTSQ(MU,X,SETA,SEPS,ALPHA,BETP,N) 15 RETURN END SUBROUTINE WLEASTSQ(MU,X,SETA,SEPS,ALPHA,BETP,N) IMPLICIT NONE REAL*8 MU(1000),X(1000),SETA,SEPS,ALPHA,BETP,W(1000) REAL*8 STAT(12),B(2),SST,SSE REAL*8 MUMATRIX(1000,2), XVECTOR(1000) INTEGER N,I,J,K C DO UNWEIGHTED REGRESSION CALL DRLINE(N,MU,X,ALPHA,BETP,STAT) C ASSIGN WEIGHTS DO 70 I=1,N W(I)=1.0/(SEPS**2+BETP**2*MU(I)**2*EXP(SETA**2)*(EXP ! (SETA**2)-1)) C STORE SQUARE ROOTS OF WEIGHTS W(I)=SQRT(W(I)) 67 CONTINUE 70 CONTINUE C CREATE MATRIX WITH WEIGHTS, WEIGHTS TIMES MU DO 99 J=1,N MUMATRIX(J,1)=W(J) MUMATRIX(J,2)=W(J)*MU(J) 99 CONTINUE C CREATE VECTOR WITH WEIGHTS TIMES X DO 101 K=1,N XVECTOR(K)=W(K)*X(K) 101 CONTINUE C FIND NEW ALPHA AND BETP CALL DRLSE(N,XVECTOR,2,MUMATRIX,1000,0,B,SST,SSE) ALPHA=B(1) BETP=B(2) RETURN END SUBROUTINE STANDDEV(X,N,SD) C THIS SUBROUTINE TAKES AN ARRAY OF VALUES AND RETURNS THE SAMPLE STANDARD DEVIATION IMPLICIT NONE REAL*8 X(1000),SD,SUM,AVG,SQSUM,SQDEV,DEV INTEGER N,I,J ! FINDING SAMPLE MEAN SUM=0 DO 10 I=1,N SUM=SUM+X(I) 10 CONTINUE AVG=SUM/N ! FINDING SUM OF SQUARES OF DEVIATIONS SQSUM=0 DO 20 J=1,N DEV=X(J)-AVG SQDEV=DEV**2 SQSUM=SQSUM+SQDEV 20 CONTINUE ! FIND STANDARD DEVIATION SD=SQRT(SQSUM/(N-1)) RETURN END SUBROUTINE DATAGROUP(MU,X,N,MU0,REP,R,MUCOUNT) IMPLICIT NONE REAL*8 MU(1000),X(1000),MU0(1000),REP(1000,1000) REAL*8 MULAST INTEGER R(1000),MUCOUNT,J,N,REPCOUNT MULAST=-1 MUCOUNT=0 DO 100 J=1,N IF (MU(J).NE.MULAST) THEN REPCOUNT=0 MUCOUNT=MUCOUNT+1 MULAST=MU(J) ENDIF REPCOUNT=REPCOUNT+1 REP(MUCOUNT,REPCOUNT)=X(J) MU0(MUCOUNT)=MU(J) R(MUCOUNT)=REPCOUNT 100 CONTINUE RETURN END SUBROUTINE SSQUARE(ALPHA,BETP,MUCOUNT,MU0,R,REP,S20) IMPLICIT NONE REAL*8 ALPHA,BETP,MU0(1000),REP(1000,1000),S20(1000) REAL*8 PRED(1000),SSQ(1000),RESID2 INTEGER R(1000),J,K,L,MUCOUNT DO 5 J=1,MUCOUNT PRED(J)=ALPHA+BETP*MU0(J) 5 CONTINUE DO 7 K=1,MUCOUNT SSQ(K)=0 DO 9 L=1,R(K) RESID2=(REP(K,L)-PRED(K))**2 SSQ(K)=SSQ(K)+RESID2 9 CONTINUE S20(K)=SSQ(K)/R(K) 7 CONTINUE RETURN END