      FUNCTION CLMST(Na,Nb,I2j,Deg)
      IMPLICIT NONE
      REAL CLMST
      INTEGER Na, Nb, I2j(2:Na+Nb)
      REAL Deg(Na+Nb)
C
C     calculates confidence level that points from sets A and B
C     are from same pdf
C
C     NA = number of points in set A
C     NB = number of points in set B
C     I2J is the array generated by a MSTREE showing the connections of
C     a minimum spanning tree constructed using the point sets A and B.
C     I2J(n) gives the number of the point to which point n is connected
C     Deg is a scratch array.  It is one less than the degree of the 
C     node, i.e., the number of edges emanating from the node, on return.
C     The points of set A are assumed to be the first NA points of the
C     combined set, followed by the points of set B.
C
c     Ref.:  Friedman, J. H. and Rafsky, L. C. "Multivariate analogs
c         of the Wald-Wolfowitz and Smirnov two-sample tests", Annals of
c         Statistics, 7, 697-717 (1979).
*
*                 W.J.Metzger    5 Jan 1999
*
      REAL ERF
      REAL c, std, SQRT2
      PARAMETER (SQRT2=1.414214)
*
*     some results of possible interest are put in common /CLMSTCM/
*
      INTEGER NINcomp
      REAL EXPct, VAR
      COMMON /CLMSTCM/ NINcomp, EXPct, VAR
*
*     printing can be controled with IPRint
*
      INTEGER IPRint
      COMMON /MSTCLFPR/ IPRint
*
      INTEGER i, n
C
C     find number of incompatible edges
C
      n = Na + Nb
      NINcomp = 0
      DO i = 2, n
         IF ( (i.LE.Na .AND. I2j(i).GT.Na) .OR. 
     &        (i.GT.Na .AND. I2j(i).LE.Na) ) NINcomp = NINcomp + 1
         Deg(i) = 0.
      ENDDO
      IF ( IPRint.NE.0 ) PRINT *, NINcomp, ' incompat edges'
C
C     find number of edges at node i (degree of node)
C
      CLMST = 0.
      DO i = 2, n
         IF ( I2j(i).GT.0 .AND. I2j(i).LE.n ) THEN
            Deg(I2j(i)) = Deg(I2j(i)) + 1.
         ELSE
            PRINT *, 'CLMST: error in I2J(', i, ' )', I2j(i)
            CLMST = -1.
         ENDIF
      ENDDO
      IF ( CLMST.LT.0. ) RETURN
C
C     Factor C for variance
C
      c = 0.
      DO i = 2, n
         c = c + Deg(i)*(Deg(i)+1.)
      ENDDO
      c = 0.5*c
C
      EXPct = FLOAT(2*Na*Nb+n)/FLOAT(n)
      VAR = FLOAT(2*Na*Nb)/FLOAT(n*n-n)
     &      *(FLOAT(2*Na*Nb-n)/n+(c-(n-2))/FLOAT((n-2)*(n-3))
     &      *(n*(n-1)-4*Na*Nb+2))
      std = ((NINcomp+1)-EXPct)/SQRT(VAR)
*
      CLMST = 0.5*(1.+ERF(std/SQRT2))
      IF ( IPRint.NE.0 ) PRINT *, 'clmst2', NINcomp + 1, Na, Nb, EXPct, 
     &                         c, std, CLMST
*
      END
      BLOCKDATA 
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      INTEGER IPRint
C*** End of declarations inserted by SPAG
      COMMON /MSTCLFPR/ IPRint
      DATA IPRint/0/
      END
