      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
