      PROGRAM NOCOM
c
c  Update history:
c   1979 originally written in Fortran IV
c   28 May 1994 ported to Microsoft Fortran
c   14 Jan 2001 ported to GNU G77 Fortran
c   16 Jan 2001 made real variables double precision
c
C     MC = maximum number of components in mixture
C     MN = maximum number of observations

C  Change the statements below wherever they occur in the program

      PARAMETER( MN = 10000 )
      PARAMETER( MC = 9 )

      implicit double precision (a-h, o-z)

      COMMON SPINIT,XM(3,MC),PROB(MC),
     1 X(MN),XX(MN),EXPOA,QJ,QQ,XLN,NS,NC,U(MC),P(MC),S(MC),
     1 FS(MC),FV(MC),UINI(MC),PINI(MC),SINI(MC),EXPO,IOPT,IU(MC),
     2 IP(MC),IS(MC),IEX,idumy,DV(MC),PS(MC),DLIM,NEWF
c  idumy required for padding. Better to have largest variables first.

      LOGICAL FOUT
C
C  IOPT=1  each component has its own standard deviation. This parametrization
C          is basically unstable, ie, the likelihood may become infinite.
C  IOPT=2  there is only one standard deviation parameter, S(1), and the
C          standard deviations of the other components are given by FS(i)*S(1),
C          i=2..NC, where the FS(i) are fixed factors.
C
      EXTERNAL F
      DATA D,TOL,RLB,RUB/0.1,0.01,-10.0,10.0/

C  DLIM = largest number such that exp(-dlim) is still nonzero
      call expmin(DLIM)

      NEWF=1

      CALL INITIAL(FOUT)

C  READ ORIG. OBS. FROM INPUT FILE
   50 i=1
   71 READ(7,*,END=72) XX(I)
      i=i+1
      if (i.le.mn) go to 71
      write(*,69) mn
   69 format(' Number of observations exceeds maximum allowed:',i8)
   72 ns=i-1
      WRITE(*,73) ns
   73 FORMAT(' Using',I8,' observations')

   79 MAL=0
   80 WRITE(*,100)
  100 FORMAT(/' ENTER EXPONENT')
      READ(*,*) EXPO
      WRITE(*,110)
  110 FORMAT(' ESTIMATE EXPONENT (1) or keep it constant (0)?')
      READ(*,*) IEX
      IF((EXPO.EQ.1.0).AND.(IEX.EQ.0)) GO TO 140
      MAL=MAL+1
      IF(MAL.GT.1) GO TO 140
      QQ=0.
      DO 130 I=1,NS
      R=XX(I)
      IF(R.GT.0.0) GO TO 130
      WRITE(*,120) I
  120 FORMAT(' OBSERVATION NO.',I8,' IS ZERO OR NEGATIVE'/)
      GO TO 999
  130 QQ=QQ+LOG(R)
  140 IF(NEWF.EQ.3) GO TO 225
      EXPOA=-100.
      WRITE(*,150) mc
  150 FORMAT(' ENTER NUMBER OF COMPONENTS (0 < n <',i2,')')
      READ(*,*) NC
      IF(NC.EQ.1) GO TO 225
      IF(NC.LE.MC) GO TO 180
      WRITE(*,160)
  160 FORMAT(' Number of components too large.'/)
      GO TO 999
  180 WRITE(*,182)
  182 FORMAT(' Fixed ratios of standard deviations (2; recommended)'
     1 /,' or separate standard deviation for each component (1)?')
      READ(*,*) IOPT
      WRITE(*,185)
  185 FORMAT(/' ENTER STARTING PARAMETER VALUES FOR EACH COMPONENT')
C
      DO 205 I=1,NC
      WRITE(*,190) I
  190 FORMAT(/' Mean for component',I2,':')
      READ(*,*) UINI(I)
      U(I)=UINI(I)
      WRITE(*,195)
  195 FORMAT(' Estimate this mean (1) or keep its value constant (0)?')
      READ(*,*) IU(I)
      IF (IOPT.NE.1) GO TO 197

c       iopt=1 -- separate standard deviations for each component
      WRITE(*,196) i
  196 FORMAT(' Standard deviation for component',I2,':')
      READ(*,*) SINI(I)
      S(I)=SINI(I)
      WRITE(*,201)
  201 FORMAT(' Estimate this std.dev. (1) or keep value constant (0)?')
      READ(*,*) IS(I)
      GO TO 199

c       iopt=2 -- fixed ratios of standard deviations
  197 if (i.eq.1) then
       WRITE(*,*) ' Standard deviation for component 1:'
       read(*,*) SINI(1)
       S(1)=SINI(1)
       fs(1)=1.0
      else
       write(*,*) ' Ratio of standard deviation to that of comp. 1:'
       READ(*,*) FS(I)
      endif

  199 WRITE(*,200) I
  200 FORMAT(' Weight (probability) of component',I2,':')
      READ(*,*) PINI(I)
      P(I)=PINI(I)
      WRITE(*,202)
  202 FORMAT(' Estimate weight (1) or keep value constant (0)?')
  205 READ(*,*) IP(I)
      SPINIT=0.0
      DO 206 I=1,NC
      IF (IP(I).EQ.0) SPINIT=SPINIT+PINI(I)
  206 CONTINUE
      SPINIT=1.0-SPINIT
      IF (SPINIT.GT.0.0) GO TO 208
      WRITE(*,207)
  207 FORMAT(' ERROR: sum of variable proportions is less than zero.'/)
      GO TO 999
  208 SPINIT=1.0/SPINIT
C
      IF (IOPT.EQ.1) GO TO 225
      WRITE(*,215)
  215 FORMAT(' Estimate standard deviations (1) or keep constant (0)?')
      READ(*,*) IS(1)
      DO 220 I=2,NC
      IS(I)=IS(1)
      FV(I)=1.0/(FS(I)*FS(I))
      SINI(I)=FS(I)*SINI(1)
  220 S(I)=SINI(I)

  225 IF(IEX.NE.0) GO TO 277

      FMAX=F(EXPO,FOUT)

      GO TO 650
  277 DD=D
      EXPO2=EXPO

      CALL FMAX1(F,EXPO2,DD,TOL,RLB,RUB,FMAX,IER2,FOUT)

      EXPO=EXPO2
      IF(IER2.EQ.0) GO TO 350
      WRITE(*,300)
  300 FORMAT(' EXPONENT OUT OF BOUNDS'/)
      GO TO 999
  350 WRITE(*,400) EXPO,FMAX
      IF(FOUT) WRITE(8,400) EXPO,FMAX
c     Output final values
  400 FORMAT(' Final values:'/' EXPONENT',F17.4/' LN(LIK)',F18.6)
C
  650 WRITE(*,660)
  660 FORMAT(/' ENTER  0 TO STOP,'/
     1 8X,'1 Start over (data will be read again)'/
     1 8X,'2 Start over, retaining current (transformed) observations'/
     1 8X,'3 New exponent only, same data, with parameter'/
     1 12X,'adjustment due to change of exponent')
      READ(*,*) NEWF
      if( (newf .lt. 0) .or. (newf .gt. 3) ) go to 650
      IF(NEWF.EQ.0) GO TO 999
      IF(NEWF.EQ.1) REWIND 7
      IF(NEWF.GE.2) GO TO 80
      GO TO 50
  999 END



      subroutine expmin(rrmax)
c  Finds largest number, rrmax, for which EXP(-rrmax) is still nonzero.
c  Assumes that EXP( -(r+1) ) is zero and does not produce an error.
      implicit double precision (a-h, o-z)
      rrmax = 1.0
  120 rrmax = rrmax + 1.0
      if (exp(-rrmax) .eq. 0.0) go to 130
      go to 120
  130 rrmax = rrmax - 1.0
      return
      end



      SUBROUTINE NOCO(FOUT)
C
C TO MAXIMIZE XLN=LOG LIKELIHOOD OVER THE PARAMETERS U=MEANS, S=STANDARD
C DEVIATION, P=PROPORTIONS IN A MIXTURE OF NORMAL DISTRIBUTIONS. IF ON INPUT,
C IER=1, THEN THE LIKELIHOOD WILL BE COMPUTED FOR THE STARTING VALUES OF THE
C PARAMETERS WITHOUT MAXIMIZATION.
C
C TT=CRITICAL DIFFERENCE IN SUCCESSIVE  LOG LIKELIHOODS.
C IEND=MAX. NO. OF ITERATIONS. IF EXCEEDED, RETURN WITH IER=1.
C NC=NUMBER OF COMPONENTS
C NS=SAMPLE SIZE
C XM(k,j) refers to the observations probabilistically assigned to the
C         j-th component, where
C XM(1,j) = S(x),  XM(2,j) = n,  XM(3,j) = S(x-xbar)**2.

      PARAMETER( MN = 10000 )
      PARAMETER( MC = 9 )

      implicit double precision (a-h, o-z)
      LOGICAL FOUT

      COMMON SPINIT,XM(3,MC),PROB(MC),
     1 X(MN),XX(MN),EXPOA,QJ,QQ,XLN,NS,NC,U(MC),P(MC),S(MC),
     1 FS(MC),FV(MC),UINI(MC),PINI(MC),SINI(MC),EXPO,IOPT,IU(MC),
     2 IP(MC),IS(MC),IEX,idumy,DV(MC),PS(MC),DLIM,NEWF

      DATA TT,IEND/.0001,300/
      XL1=-1.0E30
      ITT=0
      XS=NS
      IF(NC.GT.1) GO TO 100
C
C NC=1 (ONE COMPONENT ONLY)
      U(1)=0.
      S(1)=0.
      DO 50 I=1,NS
      U(1)=U(1)+X(I)
   50 S(1)=S(1)+X(I)*X(I)
      U(1)=U(1)/XS
      S(1)=S(1)/XS-U(1)*U(1)
      P(1)=1.
      S(1)=SQRT(S(1))
      T1=0.

      CALL LNKL

      GO TO 325
C
C NC>1. START ITERATION LOOP
  100 IF (ITT.NE.0) GO TO 274
C
C ITT=0, STARTING VALUES
      WRITE(*,271)
      IF(FOUT) WRITE(8,271)
  271 FORMAT(' STARTING VALUES')
      GO TO 300
C
C ITT>0. COMPUTE NEW PARAMETER ESTIMATES
  274 DO 275 J=1,NC
      U(J)=XM(1,J)/XM(2,J)
  275 P(J)=XM(2,J)/XS
      IF (IOPT.NE.1) GO TO 280
C IOPT=1
      DO 276 J=1,NC
  276 S(J)=SQRT(XM(3,J)/XM(2,J))
      GO TO 285
C IOPT=2
  280 S(1)=XM(3,1)
      DO 281 J=2,NC
  281 S(1)=S(1)+XM(3,J)*FV(J)
      S(1)=SQRT(S(1)/XS)
      DO 282 J=2,NC
  282 S(J)=FS(J)*S(1)
  285 DO 286 I=1,3
      DO 286 J=1,NC
      XM(I,J)=0.
  286 CONTINUE
C
C SET PARAMETERS BACK TO INITIAL VALUES WHERE DESIRED
      DO 290 J=1,NC
      IF(IU(J).EQ.0) U(J)=UINI(J)
      IF(IS(J).EQ.0) S(J)=SINI(J)
      IF(IP(J).NE.0) GO TO 288
C IP(J)=0
      P(J)=PINI(J)
      GO TO 290
C IP(J).NE.0
  288 P(J)=P(J)*SPINIT
  290 CONTINUE
C
C COMPUTE LIKELIHOOD

  300 CALL LNKL

      T1=ABS(XL1-XLN)
      IF( (MOD(ITT,25).EQ.0).OR.(T1.LE.TT) ) GO TO 325
  303 XL1=XLN
      ITT=ITT+1
      GO TO 100
C
C OUTPUT EVERY 25TH ITERATION
  325 WRITE(*,1) ITT,(U(J),J=1,NC)
      IF(FOUT) WRITE(8,1) ITT,(U(J),J=1,NC)
    1 FORMAT(1X,I3,2X,'MEANS',9F20.6)
      WRITE(*,2) (S(J),J=1,NC)
      IF(FOUT) WRITE(8,2) (S(J),J=1,NC)
    2 FORMAT(4X,'STD.DEV',9F20.6)
      WRITE(*,3) (P(J),J=1,NC)
      IF(FOUT) WRITE(8,3) (P(J),J=1,NC)
    3 FORMAT(6X,'PROP.',9F20.6)
      WRITE(*,4) EXPO,XLN
      IF(FOUT) WRITE(8,4) EXPO,XLN
    4 FORMAT(' EXPO LN(L)',F20.4,F24.6/)
C
      IF(T1-TT) 600,600,500
  500 IF(ITT-IEND) 303,510,510
  510 IER=1
  600 RETURN
      END
C  ----------


      SUBROUTINE LNKL
C  COMPUTES LOG LIKELIHOOD, XLN. FIRST, THE ABSOLUTELY SMALLEST NORMAL
C  DEVIATE IS DETERMINED (DMI). THEN, DMI IS FACTORED OUT OF THE EXPONENT.

      PARAMETER( MN = 10000 )
      PARAMETER( MC = 9 )

      implicit double precision (a-h, o-z)

      COMMON SPINIT,XM(3,MC),PROB(MC),
     1 X(MN),XX(MN),EXPOA,QJ,QQ,XLN,NS,NC,U(MC),P(MC),S(MC),
     1 FS(MC),FV(MC),UINI(MC),PINI(MC),SINI(MC),EXPO,IOPT,IU(MC),
     2 IP(MC),IS(MC),IEX,idumy,DV(MC),PS(MC),DLIM,NEWF

C INITIALIZE XLN WITH THE LOG OF THE JACOBIAN, QJ
      XLN=QJ
      DO 100 J=1,NC
  100 PS(J)=-LOG(P(J)/S(J))
C
      DO 270 K=1,NS
      DMI=1.E+30
      DO 210 J=1,NC
      D=(X(K)-U(J))/S(J)
      D=0.5*D*D+PS(J)
      IF(D.GE.DMI) GO TO 210
      DMI=D
      JJ=J
  210 DV(J)=D
      PTOT=0.
      DO 250 J=1,NC
      IF(J-JJ) 220,230,220
  220 D=DMI-DV(J)
      IF(D+DLIM) 222,225,225
  222 PROB(J)=0.
      GO TO 250
  225 PROB(J)=EXP(D)
      GO TO 240
  230 PROB(J)=1.0
  240 PTOT=PTOT+PROB(J)
  250 CONTINUE
      XLN=XLN+LOG(PTOT)-DMI
C
C COMPUTE SUMS NEEDED FOR NEW PARAMETER VALUES
      DO 270 J=1,NC
      XP=PROB(J)/PTOT
      XM(1,J)=XM(1,J)+XP*X(K)
      XM(2,J)=XM(2,J)+XP
      D=X(K)-U(J)
  270 XM(3,J)=XM(3,J)+XP*D*D
      RETURN
      END
C  ----------


      SUBROUTINE ADJM
C  ADJUSTS MEANS AND STANDARD DEVIATIONS FOR POWER TRANSFORMATION.
C  RR = OLD VALUE OF EXPONENT, SS = NEW VALUE.

      PARAMETER( MN = 10000 )
      PARAMETER( MC = 9 )

      implicit double precision (a-h, o-z)

      COMMON SPINIT,XM(3,MC),PROB(MC),
     1 X(MN),XX(MN),RR,QJ,QQ,XLN,NS,NC,U(MC),P(MC),SD(MC),
     1 FS(MC),FV(MC),UINI(MC),PINI(MC),SINI(MC),SS,IOPT

      IF(NC.EQ.1) RETURN
      R=RR
      S=SS
      IF(ABS(RR).LE.0.001) R=0.
      IF(ABS(SS).LE.0.001) S=0.
      IF(R.EQ.S) RETURN
      XB=0.
      DO 100 I=1,NC
  100 XB=XB+P(I)*U(I)
      IF(R.EQ.0.) GO TO 200
      IF(S.EQ.0.) GO TO 300
C R>0, S>0
      Q=S/R
      DO 110 I=1,NC
      SD(I)=SD(I)*((XB-R)*R+1.)**(Q-1.)
  110 U(I)=(((U(I)-R)*R+1.)**Q-1.)/S+S
      RETURN
C R=0
  200 DO 210 I=1,NC
      SD(I)=SD(I)*EXP(XB*S)
  210 U(I)=(EXP(S*U(I))-1.)/S+S
      RETURN
C S=0
  300 DO 310 I=1,NC
      SD(I)=SD(I)/((XB-R)*R+1.)
  310 U(I)=LOG((U(I)-R)*R+1.)/R
      RETURN
      END
C   ----------


      FUNCTION F(EXPO,FOUT)

      PARAMETER( MN = 10000 )
      PARAMETER( MC = 9 )

      implicit double precision (a-h, o-z)

      COMMON SPINIT,XM(3,MC),PROB(MC),
     1 X(MN),XX(MN),EXPOA,QJ,QQ,XLN,NS,NC,U(MC),P(MC),S(MC),
     1 FS(MC),FV(MC),UINI(MC),PINI(MC),SINI(MC),DUM

      LOGICAL FOUT
      DUM=EXPO
      IF(EXPO.EQ.EXPOA) GO TO 300
      IF(EXPOA.EQ.-100.) GO TO 281

      CALL ADJM

  281 CONTINUE
      QJ=0.
      IF(EXPO.EQ.1.) GO TO 290
      QJ=QQ*(EXPO-1.)
      IF(ABS(EXPO).LE.0.001) GO TO 284
      DO 282 I=1,NS
  282 X(I)=(XX(I)**EXPO-1.)/EXPO+EXPO
      GO TO 300
  284 DO 286 I=1,NS
  286 X(I)=LOG(XX(I))
      GO TO 300
  290 DO 292 I=1,NS
  292 X(I)=XX(I)
  300 CONTINUE

      CALL NOCO(FOUT)

      F=XLN
      EXPOA=EXPO
      RETURN
      END
C   ----------


      SUBROUTINE FMAX1(F,X,DEL,TOL,LB,UB,FMAX,IER,FOUT)
C  TO MAXIMIZE THE FUNCTION F OF X. WRITTEN BY R.W. COTTINGHAM.
C  F     - EXTERNAL FUNCTION OF X. MUST BE CALLABLE AS F(X).
C  X     - PARAMETER
C  DEL   - INITIAL STEP SIZE
C  TOL   - TOLERANCE OF FINAL RESULT (STOPPING CRITERION).
C  LB    - LOWER BOUND OF X
C  UB    - UPPER BOUND OF X
C  FMAX  - MAX. OF F ON OUTPUT
C  IER   - ERROR FLAG. 0=OK, 1=X EXCEEDED BOUNDS.
C  FOLD  - MAX. OF F PREVIOUS TO FMAX
C  XNEW  - NEXT VALUE OF X TRIED
C  FNEW  - VALUE OF F AT XNEW
C  NSTEP - NO. OF STEPS TAKEN IN PRESENT DIRECTION
C
      implicit double precision (a-h, o-z)
      EXTERNAL F
      double precision LB
      LOGICAL FOUT
C  INITIALIZATION
      D=DEL
      IER=1
      NSTEP=2
      FOLD=F(X,FOUT)
      X=X+D
      FMAX=F(X,FOUT)
      IF(FMAX.GE.FOLD) GO TO 10
      D=-D
      X=X+D
      TMP=FMAX
      FMAX=FOLD
      FOLD=TMP
C  SEARCH FOR MAXIMUM
   10 NSTEP=NSTEP+1
      IF(ABS(D).LE.TOL) RETURN
      XNEW=X+D
C  CHECK BOUNDARY CONDITIONS
      IF(XNEW.LT.LB .OR. XNEW.GT.UB) RETURN
      FNEW=F(XNEW,FOUT)
      IF(FNEW.LE.FMAX) GO TO 20
C        FNEW>FMAX
      FOLD=FMAX
      FMAX=FNEW
      X=XNEW
      IF(IER.EQ.1) GO TO 10
      D=-D/2.
      NSTEP=0
      GO TO 10
C        FNEW<= FMAX
   20 IER=0
      IF(NSTEP.GT.1) GO TO 30
C        NO. OF STEPS = 1
      D=-D
      IF(FNEW.GT.FOLD) FOLD=FNEW
      GO TO 10
C        NO. OF STEPS > 1
   30 D=-D/2.
      NSTEP=0
      IF(FNEW.LE.FOLD) GO TO 10
      D=-D
      FOLD=FNEW
      GO TO 10
      END



      SUBROUTINE INITIAL(FOUT)

      PARAMETER( MN = 10000 )
      PARAMETER( MC = 9 )

      implicit double precision (a-h, o-z)
      CHARACTER*1 ANSWER
      LOGICAL FOUT

C Line 1
      WRITE(*,100) CHAR(201),(CHAR(205),I=1,39),CHAR(187)
  100 FORMAT(/4X,56A1)
C Line 2
      WRITE(*,110) CHAR(186),CHAR(186)
  110 FORMAT(4X,A1,39X,A1)
C Line 3
      WRITE(*,120) CHAR(186),CHAR(186)
  120 FORMAT(4X,A1,
     1 '   NOCOM program                       ',A1)
C Line 4
      WRITE(*,130) CHAR(186),CHAR(186)
  130 FORMAT(4X,A1,
     1 '   Copyright (c) 1987-2001 Jurg Ott    ',A1)
C Line 5
      WRITE(*,110) CHAR(186),CHAR(186)
C Line 6
      WRITE(*,140) CHAR(200),(CHAR(205),I=1,39),CHAR(188)
  140 FORMAT(4X,41A1)

      WRITE(*,141) MN,MC
  141 FORMAT(/' The following maximum values are in effect:',
     1 /'   MN =',I6,' observations',
     1 /'   MC =',I6,' components')

      write(*,142)
  142 format(/' See program manual at',
     1 /'   http://linkage.rockefeller.edu/ott/nocom.htm',
     2 //' Reference the following paper when using Nocom:',
     3 /'   Ott J (1979) Detection of rare major genes in lipid levels.',
     4 /'   Hum Genet 51, 79-91'/)

      write(*,*) 'Opening "nocom.dat" file for input'
      OPEN(7,FILE='nocom.dat',STATUS='OLD',err=158)
      go to 161
  158 write(*,159)
  159 format(' File not found')
      stop

  161 write(*,162)
  162 format(/' Opening "nocom.out" file for output')
      OPEN(8,FILE='nocom.out',status='new',err=163)
      go to 168
  163 write(*,*) ' File exists. Enter Y to overwrite it,'
      write(*,*) '  N to skip writing output file, or'
      write(*,*) '  any other letter to abort program'
      read(*,*) answer
      if ( (answer.eq.'y') .or. (answer.eq.'Y') ) then
       open(8,file='nocom.out')
       fout = .true.
      else
       if ( (answer.eq.'n') .or. (answer.eq.'N') ) then
        fout = .false.
        write(*,*) '  (no output to disk file)'
       else
        stop
       end if
      end if
  168 fout = .true.

  169 IF (FOUT) WRITE(8,170)
  170 FORMAT(/6X,'Program NOCOM'/5X,
     2 ' Copyright (c) Jurg Ott 1987-2001'/
     3 1X,42('-'))

      if (fout) write(8,180)
  180 format(/' Reference the following paper when using Nocom:',
     3 /'   Ott J (1979) Detection of rare major genes in lipid levels.',
     4 /'   Hum Genet 51, 79-91'/)

      RETURN
      END
