FUNCTION AVER(X,N)
C... HAM NAY TINH TRUNG BINH SO HOC CUA CHUOI
C INPUT: + X MANG DO DAI N CHUA SO LIEU QUAN TRAC
C        + N DUNG LUONG MAU
C OUTPUT: TRUNG BINH SO HOC CUA X
C
      INTEGER N, I
      REAL X(N)		! MANG CHUA SO LIEU BAN DAU
      REAL TMP
      TMP = 0.0
      DO I = 1,N
         TMP = TMP + X(I)
      END DO
      AVER = TMP/REAL(N)
      RETURN
      END
c
      FUNCTION P_LOWER(X,N,X0)
C... HAM NAY TINH UOC LUONG XAC SUAT CUA SU KIEN
C... P(X<X0)
C INPUT: + X MANG DO DAI N CHUA SO LIEU QUAN TRAC
C.       + N DUNG LUONG MAU
C        + X0 LA MOT GIA TRI NAO DO CUA X
C OUTPUT: XAC SUAT P(X<X0)
C
      INTEGER N, N_LOWER, I
      REAL X(N)		! MANG CHUA SO LIEU BAN DAU
      N_LOWER = 0
      DO I = 1,N
         IF (X(I).LT.X0) N_LOWER = N_LOWER + 1
      END DO
      P_LOWER = REAL(N_LOWER)/REAL(N)
      RETURN
      END

	FUNCTION BERN(N,K,P)
! HAM NAY TINH XAC SUAT DE TRONG N LAN THI NGHIEM SU KIEN A 
! XUAT HIEN K LAN, VOI P=P(A)
! INPUT: + N SO LAN THI NGHIEM
!        + K SO LAN XUAT HIEN SU KIEN A
!        + P XAC SUAT CUA SU KIEN A
! OUTPUT: XAC SUAT A XUAT HIEN K NAL TRONG N LAN THI NGHIEM
!
	BERN=COMBIN(N,K)*P**K*(1-P)**(N-K)
	RETURN
	END

	FUNCTION FAC(N)
! TINH N! (N GIAI THUA)
! INPUT: + N SO NGUYEN KHONG AM
! OUTPUT: + N!
!
	REAL TMP
	  IF (N.LT.0) THEN
	     WRITE(*,*)' INVALID NUMERIC INPUT IN FAC FUNCTION'
		 STOP
	  ELSE IF (N.EQ.0.OR.N.EQ.1) THEN
	     FAC=1.0
		 RETURN
	  ELSE
	     TMP=1.0
		 DO I=2,N
	        TMP=TMP*REAL(I)
		 ENDDO
		 FAC=TMP
		 RETURN
	ENDIF
	END

      FUNCTION COMBIN(N,K)
! TINH TO HOP CHAP K CUA N
! INPUT: + N >= K LA NHUNG SO NGUYEN
! OUTPUT: TO HOP CHAP K CUA N
! FUNCTION/SUBROUTINE DUOC GOI TOI: FAC(N)
!
	IF (N.LT.0.OR.K.LT.0.OR.N.LT.K) THEN
	     WRITE(*,*)' INVALID NUMERIC INPUT IN COMBIN FUNCTION'
		 STOP
      ELSE 
         COMBIN=FAC(N)/FAC(K)/FAC(N-K)
	     RETURN
	ENDIF
	END

	FUNCTION POSSION(LAMDA,K)
! HAM NAY TINH XAC SUAT SU KIEN XUAT HIEN K LAN THEO POSSION
! INPUT: + LAMDA TRUNG BINH SO LAN XUAT HIEN
!        + K SO LAN XUA HIEN SU KIEN
! OUTPUT: XAC SUAT DE SU KIEN XUAT HIEN K LAN
!
	REAL LAMDA
	POSSION=EXP(-LAMDA)*LAMDA**K/FAC(K)
	RETURN
	END

      SUBROUTINE SORT(X,N,X1,CHISO,N1)
!
!  CHUONG TRINH NAY SAP XEP CHUOI THANH CHUOI TRINH TU
!  VA CHUOI XEP HANG
!
! INPUT: + MANG X DO DAI N CHUA CHUOI SO LIEU BAN DAU
!        + N DUNG LUONG MAU
! OUTPUT: + MANG X DO DAI N CHUA CHUOI TRINH TU
!         + MANG X1 DO DAI N1 CHUA CHUOI XEP HANG
!         + MANG CHISO DO DAI N1 CHUA CHI SO CUA CHUOI XEP HANG
!         + N1 SO THANH PHAN TRONG CHUOI XEP HANG
!
	  DIMENSION X(1),X1(1),CHISO(1)
	  INTEGER N,N1,I,J,K
	  REAL TMP
! SAP XEP THANH CHUOI TRINH TU TANG DAN
	  DO I=1,N-1
	  DO J=I+1,N
	     IF (X(J).LT.X(I)) THEN
	        TMP=X(J)
			X(J)=X(I)
			X(I)=TMP
		 ENDIF
	  ENDDO
	  ENDDO
! XAC DINH CAC THANH PHAN CUA CHUOI XEP HANG
	  N1=1
	  I=1
	  X1(N1)=X(I)
	  CHISO(N1)=REAL(I)
  10  I=I+1
      IF (I.GT.N) GOTO 100
	  IF (X(I).NE.X(I-1)) THEN
	     N1=N1+1
		 X1(N1)=X(I)
	  ENDIF
      GOTO 10
 100  CONTINUE
!  TINH CHI SO CUA CAC THANH PHAN TRONG CHUOI XEP HANG
	  DO J=1,N1
         K=0
		 CHISO(J)=0.0
		 DO I=1,N
            IF (X(I).EQ.X1(J)) THEN
			   K=K+1
			   CHISO(J)=CHISO(J)+REAL(I)
			ENDIF
		 ENDDO
		 CHISO(J)=CHISO(J)/REAL(K)
	  ENDDO
      RETURN
	END

      FUNCTION PROB1(FR,XC,M,X0)
C
C HAM NAY TINH UOC LUONG XAC SUAT P(X<X0) TU HAM PHAN BO FR
C INPUT: + FR MANG CHUA TAN SUAT TICH LUY DO DAI M
C        + XC MANG DO DAI M CHUA TRI SO CUA X MA P(X<=XC(I))=FR(XC(I))
C        + M KICH THUOC CUA FR VA XC
C        + X0 GIA TRI UNG VOI XAC SUAT CAN TINH P(X<X0)
C OUTPUT: + XAC SUAT P(X<X0)
C
	DIMENSION FR(M),XC(M)
	REAL X0
	INTEGER M,J
	IF (X0.LE.XC(1)) THEN
	   PROB1=0.0
	   RETURN
      ELSE IF (X0.GT.XC(M)) THEN
	   WRITE(*,*)' INVALID NUMERIC INPUT IN PROB1 FUNCTION'
	   STOP
	ENDIF
C TIM NHOM CHUA X0
      DO J=2,M
         IF (XC(J).GE.X0) GOTO 20
	ENDDO
  20  PROB1=FR(J-1)+(X0-XC(J-1))*(FR(J)-FR(J-1))/(XC(J)-XC(J-1))
	RETURN
	END

	SUBROUTINE FREQ(X,N,XC,FR,F,NFR,M,XMIN,STEP,INTERVAL)
C
C  CHUONG TRINH NAY TINH TAN SO (NFR) VA TAN SUAT (FR)
C  CUA MANG X DO DAI N.
C INPUT:  + MANG X CHUA GIA TRI QUAN TRAC CUA BIEN NGAU NHIEN
C         + N LA DUNG LUONG MAU
C         + INTERVAL = 0 NEU LAY SO KHOANG CHIA NGAM DINH
C                    = 1 NEU TU CHON SO KHOANG CHIA 
C OUTPUT: + M LA SO KHOANG CHIA, KHI INTERVAL = 0 THI
C           M = 5*LOG10(N)
C         + XC MANG DO DAI M CHUA TRI SO TRUNG GIAN CUA CAC KHOANG
C         + NFR MANG DO DAI M CHUA TAN SO CAC NHOM
C         + FR MANG DO DAI M CHUA TAN SUAT TICH LUY CAC NHOM (HAM PHAN BO)
C         + F MANG DO DAI M CHUA MAT DO XAC SUAT CAC NHOM
C         + XMIN GIOI HAN DUOI CUA NHOM THU NHAT (<=MIN{X(I)}
C         + STEP BUOC CHIA
C
      REAL X(1),NFR(1),FR(1),XC(1),F(1)
	INTEGER N,INTERVAL,M
	REAL XMIN,XMAX,STEP
	CHARACTER ST*10
C	   
      IF (INTERVAL.EQ.0) THEN
	     M=5*LOG10(REAL(N))
	  ELSE
 10	     WRITE(*,'(A\)') ' CHO SO KHOANG: '
		 READ*, ST
		 READ(ST,*,ERR=100) M
		 GOTO 101
 100     PRINT*,' INVALID NUMERIC INPUT ! AGAIN...'
         GOTO 10
 101     CONTINUE
      ENDIF
C
C  XAC DINH MAX, MIN CUA CHUOI DE TIM GIOI HAN CAC NHOM
      XMIN=X(1)
	XMAX=X(1)
	DO I=2,N
	   IF (X(I).LT.XMIN) XMIN=X(I)
	   IF (X(I).GT.XMAX) XMAX=X(I)
	ENDDO
C
      XMIN=AINT(XMIN)
	IF (XMIN.LT.0.0) XMIN=XMIN-1.0
	XMAX=AINT(XMAX)
	IF (XMAX.GT.0.0) XMAX=AINT(XMAX)+1.0
C
      STEP=(XMAX-XMIN)/REAL(M)
C  NEU MUON CHON BUOC TRON SO THI THEM DONG LENH SAU
C
	STEP=AINT(STEP)+1
C
C  TINH XC
      XC(1)=XMIN+STEP/2.0
	DO J=2,M
	   XC(J)=XC(J-1)+STEP
	ENDDO
C  TINH TAN SO CAC KHOANG VA HAM PHAN BO, HAM MAT DO
      STEP=STEP/2.0
	DO J=1,M
	   NFR(J)=0.0
	   DO I=1,N
	      IF ((X(I).GE.XC(J)-STEP).AND.(X(I).LT.XC(J)+STEP)) THEN
	         NFR(J)=NFR(J) + 1
	      ENDIF
	   ENDDO
	   FR(J)=NFR(J)/REAL(N)
	ENDDO
	STEP=STEP*2.0
	DO J=M,1,-1
	   F(J)=FR(J)/STEP
	   TEMP=0.0
	   DO K=1,J
	      TEMP=TEMP+FR(K)
	   ENDDO
	   FR(J)=TEMP
	ENDDO
	RETURN
	END

C
      FUNCTION PROB(FR,XCC,M,STEP,X0)
C
C HAM NAY TINH UOC LUONG XAC SUAT P(X<X0) TU HAM PHAN BO FR
C INPUT: + FR MANG CHUA TAN SUAT TICH LUY DO DAI M
C        + XC MANG CHUA TRI SO GIUA CAC NHOM DO DAI M
C        + M SO NHOM
C        + STEP CU LY NHOM (BUOC)
C        + X0 GIA TRI UNG VOI XAC SUAT CAN TINH P(X<X0)
C OUTPUT: + XAC SUAT P(X<X0)
C
	DIMENSION FR(M),XC(M),XCC(M)
	REAL STEP, X0
	INTEGER M,J
	XC=XCC
	IF (X0.LE.(XC(1)-STEP/2.0)) THEN
	   PROB=0.0
	   RETURN
      ELSE IF (X0.GE.(XC(M)+STEP/2.0)) THEN
	   PROB=1.0
	   RETURN
	ENDIF
C TIM NHOM CHUA X0
      DO J=1,M
	   XC(J)=XC(J)+STEP/2.0
         IF (XC(J).GE.X0) GOTO 20
	ENDDO
  20  PROB=FR(J-1)+(X0-XC(J-1))*(FR(J)-FR(J-1))/(XC(J)-XC(J-1))
	RETURN
	END
C
	FUNCTION X_ALFA(FR,XC,M,STEP,F0)
C
C HAM NAY TINH UOC LUONG TRI SO X0 KHI CHO TRUOC P(X<X0)=F0
C INPUT: + FR MANG CHUA TAN SUAT TICH LUY DO DAI M
C        + XC MANG CHUA TRI SO GIUA CAC NHOM DO DAI M
C        + M SO NHOM
C        + STEP CU LY NHOM (BUOC)
C        + F0=P(X<X0)=F(X0)
C OUTPUT: + TRI SO X0 UNG VOI P(X<X0)=F0
C
	DIMENSION FR(M),XC(M)
	REAL STEP, F0
	INTEGER M,J
C
	IF (F0.LE.0.0.OR.F0.GT.1.0) THEN
	    WRITE(*,*)' INVALID NUMERIC INPUT IN X_ALFA FUNCTION'
	    STOP
	ENDIF
	DO J=2,M
	   XC(J)=XC(J)+STEP/2.0
	   IF(FR(J).GE.F0) GOTO 20
	ENDDO
 20   X_ALFA=XC(J-1)+(F0-FR(J-1))/(FR(J)-FR(J-1))*(XC(J)-XC(J-1))
      RETURN
	END

	SUBROUTINE FREQ_V(X,N,XC,FR,F,NFR,M,XMAX,XMIN,STEP)
C
C  CHUONG TRINH NAY TINH TAN SO (NFR) VA TAN SUAT (FR)
C  CUA MANG X DO DAI N.
C INPUT:  + MANG X CHUA GIA TRI QUAN TRAC CUA BIEN NGAU NHIEN
C         + N LA DUNG LUONG MAU
C         + XMAX,XMIN,STEP: MAX, MIN VA CU LY NHOM
C                    
C OUTPUT: + M LA SO KHOANG CHIA
C         + XC MANG DO DAI M CHUA TRI SO TRUNG GIAN CUA CAC KHOANG
C         + NFR MANG DO DAI M CHUA TAN SO CAC NHOM
C         + FR MANG DO DAI M CHUA TAN SUAT TICH LUY CAC NHOM 
C           (HAM PHAN BO)
C         + F MANG DO DAI M CHUA MAT DO XAC SUAT CAC NHOM
C
      REAL X(1),NFR(1),FR(1),XC(1),F(1)
	INTEGER N,INTERVAL,M
	REAL XMIN,XMAX,STEP
	M=INT((XMAX-XMIN)/STEP)+1
C  TINH XC
      XC(1)=XMIN+STEP/2.0
	DO J=2,M
	   XC(J)=XC(J-1)+STEP
	ENDDO
C  TINH TAN SO CAC KHOANG VA HAM PHAN BO, HAM MAT DO
      STEP=STEP/2.0
	DO J=1,M
	   NFR(J)=0.0
	   DO I=1,N
	      IF ((X(I).GE.XC(J)-STEP).AND.(X(I).LT.XC(J)+STEP)) THEN
	         NFR(J)=NFR(J) + 1
	      ENDIF
	   ENDDO
	   FR(J)=NFR(J)/REAL(N)
	ENDDO
	STEP=STEP*2.0
	DO J=M,1,-1
	   F(J)=FR(J)/STEP
	   TEMP=0.0
	   DO K=1,J
	      TEMP=TEMP+FR(K)
	   ENDDO
	   FR(J)=TEMP
	ENDDO
	RETURN
	END

	FUNCTION Q05(X,N)
!
! HAM NAY TINH TRUNG VI CUA CHUOI X
! INPUT: + X LA MANG SO LIEU DA SAP XEP THANH CHUOI TRINH TU
!        + N LA DO DAI CUA X (DUNG LUONG MAU)
! OUTPUT: + TRUNG VI CUA CHUOI
!
      DIMENSION X(N)
	INTEGER N,I,NC
	  NC=N/2
	  IF (MOD(N,2).EQ.0) THEN
		 Q05=(X(NC)+X(NC+1))/2.0
	  ELSE
	     Q05=X(NC+1)
	  ENDIF
	  RETURN
	  END

	SUBROUTINE QUANTILES(X,N,Q50,Q25,Q75)
!
! CHUONG TRINH NAY TINH CAC PHAN VI CHINH CUA CHUOI X DO DAI N
!
! INPUT: + X MANG CHUA CHUOI SO LIEU DA SAP XEP TRINH TU
!        + N DUNG LUONG MAU
! OUTPUT: + TRUNG VI Q50
!         + CAC PHAN VI DUOI VA TREN Q25, Q75
! FUNCTION/SUBROUTINE DUOC GOI TOI: Q05(X,N)
!
      DIMENSION X(N), X1(N)
	INTEGER N,I,NC1,NC2,NT,J
	REAL Q50,Q25,Q75
	Q50=Q05(X,N)
	NT=N
	NC1=NT/2
	IF (MOD(NT,2).EQ.0) THEN
	   NC2=NC1+1
	ELSE
	   NC1=NC1+1
	   NC2=NC1
	ENDIF
	DO I=1,NC1
	   X1(I)=X(I)
	ENDDO
	Q25=Q05(X1,NC1)
	J=0
	DO I=NC2,N
	   J=J+1
	   X1(J)=X(I)
	ENDDO
	Q75=Q05(X1,NC1)
	RETURN
	END

	FUNCTION AMMGOC(X,N,K)
!
! HAM NAY TINH MOMEN GOC BAC K CUA CHUOI SO LIEU
! INPUT: + X MANG SO LIEU BAN DAU DO DAI N
!        + N DUNG LUONG MAU
!        + K BAC CUA MOMEN
! OUTPUT: MOMEN GOC BAC K
!
      DIMENSION X(N)
	INTEGER N,K,I
	REAL TMP
	TMP=0.0
	DO I=1,N
	   TMP=TMP+X(I)**K
	ENDDO
	AMMGOC=TMP/REAL(N)
	RETURN
	END


	FUNCTION AMMTAM(X,N,K)
! HAM NAY TINH MOMEN TRUNG TAM BAC K TU CHUOI SO LIEU X
! INPUT: + X MANG SO LIEU BAN DAU DO DAI N
!        + N DUNG LUONG MAU
!        + K BAC CUA MOMEN
! OUTPUT: MOMEN TRUNG TAM BAC K
! FUNCTION/SUBROUTINE DUOC GOI TOI: AMMGOC(X,N,K), COMBIN(N,K)
	DIMENSION X(N)
	INTEGER N,K,I
	REAL TMP,TB
	TB=AMMGOC(X,N,1)
	TMP=0
	DO I=0,K
	   TMP=TMP+(-1)**I*COMBIN(K,I)*TB**I*AMMGOC(X,N,K-I)
	ENDDO
	AMMTAM=TMP
	RETURN
	END

	FUNCTION BMMGOC(X,N,K)
! HAM NAY TINH MOMEN GOC BANG PHUONG PHAP PHAN NHOM
! INPUT: + X MANG SO LIEU BAN DAU DO DAI N
!        + N DUNG LUONG MAU
!        + K BAC MOMEN
! OUTPUT: MOMEN GOC BAC K
! FUNCTION/SUBROUTINE DUOC GOI TOI: 
!         FREQ(X,N,XC,FR,F,NFR,M,XMIN,STEP,INTERVAL)
!
	INTEGER N,M,I
	REAL X(N), XC(N), NFR(N), FR(N),F(N)
	REAL XMIN,STEP, TMP
	CALL FREQ(X,N,XC,FR,F,NFR,M,XMIN,STEP,0)
	TMP=0.0
	DO I=1,M
	   TMP=TMP+NFR(I)*XC(I)**K
	ENDDO
	BMMGOC=TMP/REAL(N)
	RETURN
	END

	FUNCTION BMMTAM(X,N,K)
! HAM NAY TINH MOMEN TRUNG TAM BANG PHUONG PHAP PHAN NHOM
! INPUT: + + X MANG SO LIEU BAN DAU DO DAI N
!        + N DUNG LUONG MAU
!        + K BAC MOMEN
! OUTPUT: MOMEN TRUNG TAM BAC K
! FUNCTION/SUBROUTINE DUOC GOI TOI: 
!         FREQ(X,N,XC,FR,F,NFR,M,XMIN,STEP,INTERVAL)
!         
	REAL X(N), XC(N), NFR(N), FR(N),F(N)
	INTEGER N,M,I
	REAL XMIN,STEP, TMP,TBC
	CALL FREQ(X,N,XC,FR,F,NFR,M,XMIN,STEP,0)
	TBC=0.0
	DO I=1,M
	   TBC=TBC+NFR(I)*XC(I)
	ENDDO
	TBC=TBC/REAL(N)
	DO I=1,M
	   XC(I)=XC(I)-TBC
	   TMP=TMP+NFR(I)*XC(I)**K
	ENDDO
	BMMTAM=TMP/REAL(N)
	RETURN
	END

	SUBROUTINE ANOVA(XX,N,ALFA,TB,MEDIAN,TRIMEAN,TBHC,
     &                   DX,SX,IQR,MAD,SHC,A,YULKED)
! CHUONG TRINH NAY TINH CAC DAC TRUNG THONG KE DON GIAN UNG DUNG
! TRONG PHAN TICH KHAO SAT SO LIEU
!
! INPUT: + X MANG SO LIEU BAN DAU DO DAI N
!        + N DUNG LUONG MAU
!        + ALFA (%) PHAN TRAM SO LIEU BI CAT BO
! OUTPUT: + TB TRUNG BINH SO HOC
!         + MEDIAN TRUNG VI
!         + TRIMEAN
!         + TBHC TRUNG BINH HIEU CHINH
!         + DX PHUONG SAI
!         + SX DO LECH CHUAN
!         + IQR CHI SO BIEN DO PHAN TU
!         + MAD CHI SO MAD (MEDIAN CUA DO LECH TUYET DOI SO VOI
!                            MEDIAN BAN DAU CUA CHUOI)
!         + SHC PHUONG SAI HIEU CHINH
!         + A DO BAT DOI XUNG
!         + YULKED CHI SO YULE-KENDALL
! FUNCTION/SUBROUTINE DUOC GOI TOI: SORT, QUANTILES
!
	DIMENSION XX(N),X(N),X1(N),CHISO(N)
	REAL ALFA,TB,MEDIAN,TRIMEAN,TBHC
	REAL DX,SX,IQR,MAD,SHC,A,YULKED
	REAL Q50,Q25,Q75,TMP
	INTEGER N,I,N1,K
!
	X=XX
	TB=AMMGOC(X,N,1)
	DX=AMMTAM(X,N,2)
	SX=SQRT(DX)
	A =AMMTAM(X,N,3)/(DX*SX)

	CALL SORT(X,N,X1,CHISO,N1)
	CALL QUANTILES(X,N,Q50,Q25,Q75)
!
	MEDIAN=Q50
	TRIMEAN=(Q25+2.0*Q50+Q75)/4.0
	IQR=Q75-Q25
	YULKED=(Q25-2.0*Q50+Q75)/IQR
!
	K=INT(REAL(N)*ALFA)
	TMP=0.0
	DO I=K+1,N-K
	   TMP=TMP+X(I)
	ENDDO
	TBHC=TMP/REAL(N-2*K)
!
	TMP=0.0
	DO I=K+1,N-K
	   TMP=TMP+(X(I)-TBHC)**2
	ENDDO
	SHC=TMP/REAL(N-2*K)
!
	DO I=1,N
	   X(I)=ABS(X(I)-Q50)
	ENDDO
	CALL QUANTILES(X,N,Q50,Q25,Q75)
	MAD=Q50
	RETURN
	END

	SUBROUTINE STATISTICS(XX,N,M,ALFA,TB,DX,SX,A,E,MEDIAN,
     &             TRIMEAN,TBHC,IQR,MAD,SHC,YULKED)
! CHUONG TRINH NAY TINH CAC DAC TRUNG THONG	KE CUA HE CAC 
! DAI LUONG NGAU NHIEN MA X LA MA TRAN SO LIEU BAN DAU
! INPUT: + X(N,M) MA TRAN SO LIEU N HANG (N QUAN TRAC), M COT
!          (M BIEN)
!        + N DUNG LUONG MAU
!        + M BIEN
!        + ALFA (%) SO PHAN TRAM CAC THANH PHAN BI CAT BO
! OUTPUT: + TB,DX,SX,A,E LA CAC MANG DO DAI M TUONG UNG CHUA
!           TRUNG BINH, PHUONG SAI, DO LECH CHUAN, DO BAT DOI
!           XUNG VA DO NHON CUA CAC BIEN
!         + TRIMEAN,TBHC,IQR,MAD,SHC,YULKED LA CAC MANG DO DAI
!           M TUONG UNG CHUA TRIMEAN, TRUNG BINH HIEU CHINH,
!           BIEN DO PHAN TU, MAD (MEDIAN ABSOLUTE DEVIATION),
!           PHUONG SAI HIEU CHINH VA CHI SO YULE-KENDALL CUA
!           CAC BIEN
!         + MEDIAN MANG DO DAI M CHUA TRUNG VI
!
	DIMENSION X(N,M),XX(N*M)
	INTEGER N,M,I,J,K
	DIMENSION TB(M),DX(M),SX(M),A(M),E(M)
	REAL TRIMEAN(M),TBHC(M),IQR(M),MAD(M),SHC(M),YULKED(M)
	REAL MEDIAN(M)
	DIMENSION XTMP(N)
	K=0
	DO J=1,M
	  DO I=1,N
	     K=K+1
		 X(I,J)=XX(K)
	  ENDDO
	ENDDO
	DO 10 J=1,M
	   DO I=1,N
	      XTMP(I)=X(I,J)
	   ENDDO
	   TB(J)=AMMGOC(XTMP,N,1)
	   DX(J)=AMMTAM(XTMP,N,2)
	   SX(J)=SQRT(DX(J))
	   A(J)=AMMTAM(XTMP,N,3)/(DX(J)*SX(J))
	   E(J)=AMMTAM(XTMP,N,4)/(DX(J)*DX(J))-3.0
	   CALL ANOVA(XTMP,N,ALFA,TB(J),MEDIAN(J),TRIMEAN(J),TBHC(J),
     &            DX(J),SX(J),IQR(J),MAD(J),SHC(J),A(J),YULKED(J))
  10  CONTINUE
	RETURN
	END

	SUBROUTINE BINOMD(N,P,X,PX,FX)
! CHUONG TRINH NAY TINH PHAN BO NHI THUC
! INPUT: + N SO LAN THI NGHIEM
!        + P=P(A) XAC SUAT SU KIEN A
! OUTPUT: + X MANG DO DAI N+1 (TU 0..N) CHUA GIA TRI CUA X
!         + PX MANG DO DAI N+1 (TU 0..N) CHUA XAC SUAT
!           DE X NHAN CAC GIA TRI 0..N
!         + FX MANG DO DAI N+1 (TU 0..N) CHUA XAC SUAT TICH
!           LUY CUA X (FX(K)=P(X<K))
	DIMENSION X(0:N),PX(0:N),FX(0:N)
	DO K=0,N
	   X(K)=K
	   PX(K)=BERN(N,K,P)
	ENDDO
	DO K=0,N
	   FX(K)=0.0
	   DO I=0,K
	      FX(K)=FX(K)+PX(I)
	   ENDDO
      ENDDO
	RETURN
	END

	SUBROUTINE POSSD(N,LAMDA,X,PX,FX)
! CHUONG TRINH NAY TINH PHAN BO POSSION
! INPUT: + N SO LAN THI NGHIEM
!        + LAMDA TRUNG BINH SO LAN XUAT HIEN
! OUTPUT: + X MANG DO DAI N+1 (TU 0..N) CHUA GIA TRI CUA X
!         + PX MANG DO DAI N+1 (TU 0..N) CHUA XAC SUAT
!           DE X NHAN CAC GIA TRI 0..N
!         + FX MANG DO DAI N+1 (TU 0..N) CHUA XAC SUAT TICH
!           LUY CUA X (FX(K)=P(X<K))
	DIMENSION X(0:N),PX(0:N),FX(0:N)
	REAL LAMDA

	DO K=0,N
	   X(K)=K
	   PX(K)=POSSION(LAMDA,K)
	ENDDO
	DO K=0,N
	   FX(K)=0.0
	   DO I=0,K
	      FX(K)=FX(K)+PX(I)
	   ENDDO
      ENDDO
	RETURN
	END

	FUNCTION FNORMD(TB,SX,X0)
! HAM NAY TINH GIA TRI HAM PHAN BO CUA BIEN NGAU NHIEN 
! PHAN BO CHUAN CO KY VONG = TB, PHUONG SAI =SX**2
! INPUT: + TB KY VONG CUA X
!        + SX DO LECH BINH PHUONG TRUNG BINH CUA X
!        + X0 GIA TRI NAO DO CUA X
! OUTPUT: FNORMD = F(X<X0)
!
	REAL TB,SX,X,P1
	IF (TB.EQ.X0) THEN
	   P1=0.5
	ELSE IF (TB.GT.X0) THEN
	   P1=0.5-ANORMD(TB,SX,X0,TB)
	ELSE
	   P1=0.5+ANORMD(TB,SX,TB,X0)
	ENDIF
	FNORMD=P1
	RETURN
	END

	FUNCTION ANORMD(AMUY,XICMA,A,B)
! HAM NAY TINH XAC SUAT DE BIEN NGAU NHIEN PHAN BO CHUAN
! NHAN GIA TRI TRONG KHOANG (A,B): P(A<X<B). THUC CHAT
! HAM NAY LA TINH TICH PHAN TU A-B
! INPUT: + MUY KY VONG TOAN
!        + XICMA DO LECH BINH PHUONG TRUNG BINH
!        + A VA B LA HAI CAN TICH PHAN.
! OUTPUT: P(A<X<B)
!
	REAL AMUY,XICMA,A,B,C
	PARAMETER(PI=3.1415927,EP=1.0E-6)
	F(X)=1.0/(SQRT(2.0*PI)*XICMA)*EXP(-0.5*((X-AMUY)/XICMA)**2)
	IF (A.EQ.B) THEN
	   ANORMD=0.0
	   RETURN
	ENDIF
	IF (A.GT.B) THEN
	   C=A
	   A=B
	   B=C
	ENDIF
	S=0.0
	N=0
  10  N=N+1
      DELX=(B-A)/REAL(N)
	S1=0.0
	DO I=0,N-1
	   X=A+DELX*REAL(I)
	   F1=F(X)
	   X=X+DELX
	   F2=F(X)
	   S1=S1+(F1+F2)*DELX/2.0
	ENDDO
	SS=ABS((S1-S)/S1)
	IF (SS.GT.EP) THEN
	   S=S1
	   GOTO 10
	ENDIF
      ANORMD=S1
	RETURN
	END

	FUNCTION ANINV(P,MUY,XICMA)
! HAM NAY TINH GIA TRI X0 CUA BIEN NGAU NHIEN X PHAN BO CHUAN
! VOI HAI THAM SO MUY VA XICMA KHI P(X<X0)=P.
! KHI TINH SE SU DUNG BIEN PHAN BO CHUAN CHUAN HOA
! INPUT: + P XAC SUAT DE X<X0 (P(X<X0)=P)
!        + MUY KY VONG CUA X
!        + XICMA DO LECH BINH PHUONG TRUNG BINH
! OUTPUT: X0
!
	PARAMETER (EPS=1.0E-6)
	REAL P,AP,MUY,XICMA,A,B,C,P0
	INTEGER IFLAG
	IF (P.EQ.0.0) THEN
	   WRITE(*,*)' THE P VALUE TOO SMALL IN ANINV'
	   STOP
	ELSE IF (P.EQ.0.5) THEN
	   ANINV=MUY
	   RETURN
	ENDIF
!
	IF (P.LT.0.5) THEN
	   AP=0.5-P
	   IFLAG=1
	ELSE IF (P.GT.0.5) THEN
	   AP=P-0.5
	   IFLAG=0
	   PRINT*,P,A,B
	ENDIF
	A=0.0
	B=6.0
	C=B

 10	P0=ANORMD(0.0,1.0,A,B)
	SS=(ABS(P0-AP)/AP)
	IF (SS.GE.EPS) THEN
	   IF (AP.LT.P0) THEN
	      C=B
	      B=(A+B)/2.0
		  GOTO 10
	   ELSE IF (AP.GT.P0) THEN
		  B=(C+B)/2.0
		  GOTO 10
	   ENDIF
	ENDIF
	A=(C+B)/2.0
	IF (IFLAG.EQ.1) THEN
 	   ANINV=-A*XICMA+MUY
	ELSE
         ANINV=A*XICMA+MUY
	ENDIF
	RETURN
	END

	FUNCTION CHINV(P,N)
! HAM NAY TINH GIA TRI X0 CUA BIEN NGAU NHIEN X PHAN BO 
! "KHI BINH PHUONG" VOI N BAC TU DO THOA MAN DIEU KIEN 
! P(X>X0)=P.
! INPUT: + P XAC SUAT DE X>X0 (P(X>X0)=P)
!        + N SO BAC TU DO, PHAI LA MOT SO THUC
! OUTPUT: X0
! SUBROUTINE/FUNCTION DUOC GOI TOI: CHIDIST
	PARAMETER (EPS=1.0E-6)
	REAL P,AP,N, A,B,C,P0

	IF (P.LT.0.0.OR.P.GT.1.0) THEN
	   WRITE(*,*)' INVALID NUMERIC INPUT IN AKBPINV FUNCTION'
	   STOP
	ENDIF
	IF (P.EQ.0.0) THEN
	   CHINV=0.0
	   RETURN
	ELSE IF (P.EQ.1) THEN
	   WRITE(*,*)' THE P VALUE TOO BIG IN AKBPINV FUNCTION'
	   STOP
	ENDIF
!
	A=0.0
	B=9999999.0
	C=A
	AP=P

 10	P0=CHIDIST(B,N)
	SS=ABS((P0-AP)/AP)
	SS1=ABS((B-C)/B)
	IF (SS.GE.EPS.AND.SS1.GE.EPS) THEN
	   IF (AP.GT.P0) THEN
	      C=B
	      B=(A+B)/2.0
		  GOTO 10
	   ELSE IF (AP.LT.P0) THEN
	      B=(C+B)/2.0
		  GOTO 10
	   ENDIF
	ENDIF
	A=(C+B)/2.0
      CHINV=A
	RETURN
	END


	FUNCTION CHIDIST(X0,N)
! HAM NAY TINH XAC SUAT P=P(X>0) VOI X LA BIEN NGAU NHIEN
! CO PHAN BO "KHI BINH PHUONG" VOI N BAC TU DO.
! INPUT: + X0 MOT GIA TRI NAO DO CUA X
!        + N BAC TU DO
! OUTPUT: P=P(X>X0)
! SUBROUTINE/FUNCTION DUOC GOI TOI: GAMMQ
      REAL X0,N
	CHIDIST=GAMMQ(N/2.0,X0/2.0)
	RETURN
	END

      FUNCTION gammq(a,x)
      REAL a,gammq,x
!CU    USES gcf,gser
      REAL gammcf,gamser,gln
      if(x.lt.0..or.a.le.0.)pause 'bad arguments in gammq'
      if(x.lt.a+1.)then
        call gser(gamser,a,x,gln)
        gammq=1.-gamser
      else
        call gcf(gammcf,a,x,gln)
        gammq=gammcf
      endif
      return
      END

      SUBROUTINE gcf(gammcf,a,x,gln)
      INTEGER ITMAX
      REAL a,gammcf,gln,x,EPS,FPMIN
      PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30)
!CU    USES gammln
      INTEGER i
      REAL an,b,c,d,del,h,gammln
      gln=gammln(a)
      b=x+1.-a
      c=1./FPMIN
      d=1./b
      h=d
      do 11 i=1,ITMAX
        an=-i*(i-a)
        b=b+2.
        d=an*d+b
        if(abs(d).lt.FPMIN)d=FPMIN
        c=b+an/c
        if(abs(c).lt.FPMIN)c=FPMIN
        d=1./d
        del=d*c
        h=h*del
        if(abs(del-1.).lt.EPS)goto 1
11    continue
      pause 'a too large, ITMAX too small in gcf'
1     gammcf=exp(-x+a*log(x)-gln)*h
      return
      END


      SUBROUTINE gser(gamser,a,x,gln)
      INTEGER ITMAX
      REAL a,gamser,gln,x,EPS
      PARAMETER (ITMAX=100,EPS=3.e-7)
!CU    USES gammln
      INTEGER n
      REAL ap,del,sum,gammln
      gln=gammln(a)
      if(x.le.0.)then
        if(x.lt.0.)pause 'x < 0 in gser'
        gamser=0.
        return
      endif
      ap=a
      sum=1./a
      del=sum
      do 11 n=1,ITMAX
        ap=ap+1.
        del=del*x/ap
        sum=sum+del
        if(abs(del).lt.abs(sum)*EPS)goto 1
11    continue
      pause 'a too large, ITMAX too small in gser'
1     gamser=sum*exp(-x+a*log(x)-gln)
      return
      END

      FUNCTION gammln(xx)
      REAL gammln,xx
      INTEGER j
      DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
      SAVE cof,stp
      DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, 
     & 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, 
     & -.5395239384953d-5,2.5066282746310005d0/
      x=xx
      y=x
      tmp=x+5.5d0
      tmp=(x+0.5d0)*log(tmp)-tmp
      ser=1.000000000190015d0
      do 11 j=1,6
        y=y+1.d0
        ser=ser+cof(j)/y
11    continue
      gammln=tmp+log(stp*ser/x)
      return
      END

	FUNCTION TINV(P,N)
! HAM NAY TINH GIA TRI X0 CUA BIEN NGAU NHIEN X PHAN BO 
! STUDENT VOI N BAC TU DO THOA MAN DIEU KIEN 
! P(|X|>X0)=P.
! INPUT: + P XAC SUAT DE |X|>X0 (P(|X|>X0)=P)
!        + N SO BAC TU DO, PHAI LA MOT SO THUC
! OUTPUT: X0
! SUBROUTINE/FUNCTION DUOC GOI TOI: TDIST
	PARAMETER (EPS=1.0E-6)
	REAL P,AP,N, A,B,C,P0
	IF (P.LT.0.0.OR.P.GT.1.0) THEN
	   WRITE(*,*)' INVALID NUMERIC INPUT IN AKBPINV FUNCTION'
	   STOP
	ENDIF
	IF (P.EQ.0.0) THEN
	   WRITE(*,*)' THE P VALUE TOO SMALL IN AKBPINV FUNCTION'
	   STOP
	ELSE IF (P.EQ.1) THEN
	   WRITE(*,*)' THE P VALUE TOO BIG IN AKBPINV FUNCTION'
	   STOP
	ENDIF
!
	A=0.0
	B=99999.0
	AP=P

 10	P0=TDIST(B,N)
	SS=ABS((P0-AP)/AP)
	SS1=ABS((B-C)/B)
	IF (SS.GE.EPS.AND.SS1.GE.EPS) THEN
	   IF (AP.GT.P0) THEN
	      C=B
		  B=(A+B)/2.0
	      GOTO 10
	   ELSE IF (AP.LT.P0) THEN
	      B=(C+B)/2.0
	      GOTO 10
	   ENDIF
	ENDIF
	A=(C+B)/2.0
      TINV=A
	RETURN
	END


	FUNCTION TDIST(X0,N)
! HAM NAY TINH XAC SUAT DE BIEN NGAU NHIEN X CO PHAN BO STUDENT
! VOI N BAC TU DO NHAN GIA TRI NGOAI KHOANG (-X0;X0):
! P=P(|X|>X0)=1-P(|X|<X0)
! INPUT: + X0 MOT GIA TRI NAO DO CUA X
!        + N SO BAC TU DO, PHAI LA MOT SO THUC
! SUBROUTINE/FUNCTION DUOC GOI TOI: BETAI
!
	REAL X0,N,X
	X=N/(N+X0**2)
	TDIST=BETAI(N/2.0,0.5,X)
	RETURN
	END 

      FUNCTION betai(a,b,x)
      REAL betai,a,b,x
!CU    USES betacf,gammln
      REAL bt,betacf,gammln
      if(x.lt.0..or.x.gt.1.)pause 'bad argument x in betai'
      if(x.eq.0..or.x.eq.1.)then
        bt=0.
      else
        bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.-x))
      endif
      if(x.lt.(a+1.)/(a+b+2.))then
        betai=bt*betacf(a,b,x)/a
        return
      else
        betai=1.-bt*betacf(b,a,1.-x)/b
        return
      endif
      END

      FUNCTION betacf(a,b,x)
      INTEGER MAXIT
      REAL betacf,a,b,x,EPS,FPMIN
      PARAMETER (MAXIT=100,EPS=3.e-7,FPMIN=1.e-30)
      INTEGER m,m2
      REAL aa,c,d,del,h,qab,qam,qap
      qab=a+b
      qap=a+1.
      qam=a-1.
      c=1.
      d=1.-qab*x/qap
      if(abs(d).lt.FPMIN)d=FPMIN
      d=1./d
      h=d
      do 11 m=1,MAXIT
        m2=2*m
        aa=m*(b-m)*x/((qam+m2)*(a+m2))
        d=1.+aa*d
        if(abs(d).lt.FPMIN)d=FPMIN
        c=1.+aa/c
        if(abs(c).lt.FPMIN)c=FPMIN
        d=1./d
        h=h*d*c
        aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2))
        d=1.+aa*d
        if(abs(d).lt.FPMIN)d=FPMIN
        c=1.+aa/c
        if(abs(c).lt.FPMIN)c=FPMIN
        d=1./d
        del=d*c
        h=h*del
        if(abs(del-1.).lt.EPS)goto 1
11    continue
      pause 'a or b too big, or MAXIT too small in betacf'
1     betacf=h
      return
      END

	FUNCTION FINV(P,N1,N2)
! HAM NAY TINH GIA TRI X0 CUA BIEN NGAU NHIEN X PHAN BO 
! FISHER (F) VOI N1 VA N2 BAC TU DO THOA MAN DIEU KIEN 
! P(X>X0)=P.
! INPUT: + P XAC SUAT DE X>X0 (P(X>X0)=P)
!        + N1 VA N2 LA SO BAC TU DO, PHAI LA NHUNG SO THUC
! OUTPUT: X0
! SUBROUTINE/FUNCTION DUOC GOI TOI: FDIST
	PARAMETER (EPS=1.0E-6)
	REAL P,AP,N1,N2, A,B,C,P0,SS,SS1
	IF (P.LT.0.0.OR.P.GT.1.0) THEN
	   WRITE(*,*)' INVALID NUMERIC INPUT IN FINV FUNCTION'
	   STOP
	ENDIF
	IF (P.EQ.0.0) THEN
	   WRITE(*,*)' THE P VALUE TOO SMALL IN FINV FUNCTION'
	   STOP
	ELSE IF (P.EQ.1) THEN
	   FINV=0.0
	   RETURN
	ENDIF
!
	A=0.0
	B=9999999.0
	C=A
	AP=P

 10	P0=FDIST(B,N1,N2)
	SS=ABS((P0-AP)/AP)
	SS1=ABS((B-C)/B)
	IF (SS.GE.EPS.AND.SS1.GE.EPS) THEN
	   IF (AP.GT.P0) THEN
	      C=B
		  B=(A+B)/2.0
		  GOTO 10
	   ELSE IF (AP.LT.P0) THEN
		  B=(C+B)/2.0
		  GOTO 10
	   ENDIF
	ENDIF
	A=(C+B)/2.0
      FINV=A
	RETURN
	END


	FUNCTION FDIST(X0,N1,N2)
! HAM NAY TINH XAC SUAT DE BIEN NGAU NHIEN X CO PHAN BO FISHER
! VOI N1 VA N2 BAC TU DO NHAN GIA TRI >X0:
! P=P(X>X0)=1-P(X<X0)
! INPUT: + X0 MOT GIA TRI NAO DO CUA X
!        + N1,N2 SO BAC TU DO, PHAI LA NHUNG SO THUC
! SUBROUTINE/FUNCTION DUOC GOI TOI: BETAI
!
	REAL X0,N1,N2,X
	X=N2/(N2+N1*X0)
	FDIST=BETAI(0.5*N2,0.5*N1,X)
	RETURN
	END 

	SUBROUTINE CORRE(X,N,M,R)
! CHUONG TRINH NAY TINH MA TRAN TUONG QUAN CHUAN HOA
! CUA TAP M BIEN TU SO LIEU BAN DAU CO DUNG LUONG N
! INPUT: + X MANG MOT CHIEU KICH THUOC N*M LUU TRU SO
!          LIEU BAN DAU DANG MA TRAN N HANG M COT
!          (X(1,1), X(2,1),...,X(N,1),X(1,2),...)
!        + N DUNG LUONG MAU
!        + M SO BIEN
! OUTPUT: + R MANG MOT CHIEU KICH THUOC M*M LUU TRU MA TRAN
!            TUONG QUAN CHUAN HOA CUA M BIEN
!
	DIMENSION X(N*M),XX(N,M),R(M*M),TB(M),SX(M)
	K=0
	DO J=1,M
	   TB(J)=0.0
	   SX(J)=0.0
	   DO I=1,N
	      K=K+1
	      TB(J)=TB(J)+X(K)
	      SX(J)=SX(J)+X(K)*X(K)
	      XX(I,J)=X(K)
	   ENDDO
	   TB(J)=TB(J)/REAL(N)
	   SX(J)=SQRT(SX(J)/REAL(N)-TB(J)*TB(J))
	ENDDO
	JK=0
	DO J=1,M
	   DO K=1,M
	      JK=JK+1
	      R(JK)=0.0
	      DO I=1,N
	         R(JK)=R(JK)+XX(I,J)*XX(I,K)
	      ENDDO
	      R(JK)=R(JK)/REAL(N)-TB(J)*TB(K)
	      R(JK)=R(JK)/(SX(J)*SX(K))
	   ENDDO
	ENDDO
	RETURN
	END

	SUBROUTINE ULKHOANGR(R,ALFA,N,R1,R2)
! CHUONG TRINH NAY TINH UOC LUONG KHOANG CUA HE SO TUONG QUAN
! (R1, R2) ==> R1<RO<R2
! INPUT: + R HE SO TUONG QUAN MAU
!        + ALFA MUC Y NGHIA
!        + N DUNG LUONG MAU
! OUTPUT: + R1 GIOI HAN DUOI CUA KHOANG
!         + R2 GIOI HAN TREN CUA KHOANG
!
	UALFA=ANINV(1.0-ALFA/2.0,0.0,1.0)
	Z=0.5*LOG((1.0+R)/(1.0-R))
	A=0.5*R/(REAL(N-1))
	B=UALFA/SQRT(REAL(N-3))
	R1=Z-A-B
	R2=Z-A+B
!	PRINT*,Z,A,B,R1,R2
	A=EXP(2.0*R1)
	R1=(A-1.0)/(A+1.0)
	A=EXP(2.0*R2)
	R2=(A-1.0)/(A+1.0)
	RETURN
	END

	SUBROUTINE COVAR(X,N,M,R)
! CHUONG TRINH NAY TINH MA TRAN TUONG QUAN CUA M BIEN
! INPUT: + X MANG MOT CHIEU KICH THUOC N*M CHUA SO LIEU BAN DAU
!          LUU THEO COT CUA MA TRAN XX(N,M)
!        + N DUNG LUONG MAU
!        + M SO BIEN
! OUTPUT: R MANG MOT CHIEU KICH THUOC M*M CHUA MA TRAN TUONG QUAN
!
	DIMENSION X(N*M),XX(N,M),R(M*M),RR(M,M),TB(M)

	K=0
	DO J=1,M
	   TB(J)=0.0
	   DO I=1,N
	      K=K+1
	      TB(J)=TB(J)+X(K)
	      XX(I,J)=X(K)
	   ENDDO
	   TB(J)=TB(J)/REAL(N)
	ENDDO
	DO J=1,M
	   DO K=J,M
	      RR(J,K)=0.0
	      DO I=1,N
	         RR(J,K)=RR(J,K)+(XX(I,J)-TB(J))*(XX(I,K)-TB(K))
	      ENDDO
	      RR(J,K)=RR(J,K)/REAL(N)
	      RR(K,J)=RR(J,K)
	   ENDDO
	ENDDO
	L=0
	DO K=1,M
	   DO J=1,M
	      L=L+1
		  R(L)=RR(J,K)
	   ENDDO
	ENDDO	
	RETURN
	END

	FUNCTION TQHANG(X,Y,N)
	DIMENSION X(N),Y(N),XX(N),YY(N),U(N),V(N)
	XX=X
	YY=Y
	DO I=1,N-1
	   DO J=I+1,N
	      IF (XX(J).LT.XX(I)) THEN
	         TMP=XX(I)
	         XX(I)=XX(J)
	         XX(J)=TMP
	      ENDIF
	      IF (YY(J).LT.YY(I)) THEN
	         TMP=YY(I)
	         YY(I)=YY(J)
	         YY(J)=TMP
	      ENDIF
	   ENDDO
	ENDDO
	DO I=1,N
	   J=0
   10    J=J+1
	   IF (X(I).EQ.XX(J)) THEN
	      U(I)=J
	      GOTO 11
	   ENDIF
	   IF (J.LT.N) GOTO 10
   11    CONTINUE
	ENDDO

	DO I=1,N
	   J=0
   12    J=J+1
	   IF (Y(I).EQ.YY(J)) THEN
	      V(I)=J
	      GOTO 13
	   ENDIF
	   IF (J.LT.N) GOTO 12
   13    CONTINUE
	ENDDO
	TMP=0.0
	DO I=1,N
!	WRITE(*,'(2F4.0)') U(I),V(I)
	   TMP=TMP+(U(I)-V(I))**2
	ENDDO
	TQHANG=1.0-6.0*TMP/REAL(N*(N-1)*(N+1))
	RETURN
	END

	SUBROUTINE HQTT_1(X,Y,N,A0,A1,R,Q,U,S,F)
	DIMENSION X(N),Y(N)
	TBX=0.0
	TBY=0.0
	DO I=1,N
	   TBX=TBX+X(I)
	   TBY=TBY+Y(I)
	ENDDO
	TBX=TBX/REAL(N)
	TBY=TBY/REAL(N)
	XY=0.0
	XX=0.0
	YY=0.0
	DO I=1,N
	   TX=X(I)-TBX
	   TY=Y(I)-TBY
	   XY=XY+TX*TY
	   XX=XX+TX*TX
	   YY=YY+TY*TY
	ENDDO
	A1=XY/XX
	A0=TBY-A1*TBX
	U=A1*XY
	Q=YY-U
	R=SQRT(U/YY)
	S=SQRT(Q/REAL(N-2))
	F=U*REAL(N-2)/Q
	RETURN
	END

	FUNCTION HSTQ(X,Y,N)
! HAM NAY TINH HE SO TUONG QUAN GIUA HAI BIEN X VA Y
! INPUT: + X MANG MOT CHIEU KICH THUOC N CHUA SO LIEU CUA X
!        + Y MANG MOT CHIEU KICH THUOC N CHUA SO LIEU CUA Y
!        + N DUNG LUONG MAU
! OUTPUT: HE SO TUONG QUAN R GIUA X VA Y
!
	DIMENSION X(N),Y(N)
	TBX=0.0
	TBY=0.0
	SX=0.0
	SY=0.0
	RXY=0.0
	DO I=1,N
	   TBX=TBX+X(I)
	   TBY=TBY+Y(I)
	   SX=SX+X(I)*X(I)
	   SY=SY+Y(I)*Y(I)
	   RXY=RXY+X(I)*Y(I)
	ENDDO
	TBX=TBX/REAL(N)
	TBY=TBY/REAL(N)
	SX=SQRT(SX/REAL(N))
	SY=SQRT(SY/REAL(N))
	RXY=RXY/REAL(N)-TBX*TBY
	HSTQ=RXY/SX/SY
	RETURN
	END

	FUNCTION TYSOTQ(X,Y,N)
! CHUONG TRINH NAY TINH TY SO TUONG QUAN GIUA HAI BIEN X VA Y
! INPUT: + X MANG MOT CHIEU KICH THUOC N CHUA SO LIEU CUA X
!        + Y MANG MOT CHIEU KICH THUOC N CHUA SO LIEU CUA Y
!        + N DUNG LUONG MAU
! OUTPUT: TY SO TUONG QUAN ETA GIUA X VA Y
!
	DIMENSION X(N),Y(N),TSO(N/2,N/2)
	DIMENSION XC(N/2),YC(N/2),TSX(N/2),TSY(N/2),A(N/2)
	LOGICAL LX,LY

	XMIN=X(1)
	XMAX=X(1)
	YMIN=Y(1)
	YMAX=Y(1)
	DO I=2,N
	   IF (X(I).LT.XMIN) XMIN=X(I)
	   IF (X(I).GT.XMAX) XMAX=X(I)
	   IF (Y(I).LT.YMIN) YMIN=Y(I)
	   IF (Y(I).GT.YMAX) YMAX=Y(I)
	ENDDO
	XMIN=1.0*INT(XMIN)
	YMIN=1.0*INT(YMIN)
	XMAX=1.0*INT(XMAX)+1.0
	YMAX=1.0*INT(YMAX)+1.0

	NSTEP=INT(5.0*LOG10(REAL(N)))+1
	STEPX=(XMAX-XMIN)/REAL(NSTEP)
	STEPY=(YMAX-YMIN)/REAL(NSTEP)

! TINH BANG PHAN BO TAN SO HAI CHIEU (TSO)
	DO KX=1,NSTEP
	   AX=XMIN+REAL(KX-1)*STEPX
	   BX=AX+STEPX

	   DO KY=1,NSTEP
	      AY=YMIN+REAL(KY-1)*STEPY
	      BY=AY+STEPY
	      TSO(KX,KY)=0.0

	      DO I=1,N
	         LX=(X(I).GE.AX).AND.(X(I).LT.BX)
	         LY=(Y(I).GE.AY).AND.(Y(I).LT.BY)
	         IF (LX.AND.LY) THEN
	            TSO(KX,KY)=TSO(KX,KY)+1.0
	         ENDIF
	      ENDDO
         ENDDO
	ENDDO

	XMIN=XMIN-STEPX/2.0
	YMIN=YMIN-STEPY/2.0

! TINH CAC TRI SO GIUA CUA X (XC) VA Y (YC)
	DO KX=1,NSTEP
	   XC(KX)=XMIN+REAL(KX)*STEPX
	   TSX(KX)=0.0
	   DO KY=1,NSTEP
	      TSX(KX)=TSX(KX)+TSO(KX,KY)
	   ENDDO
	ENDDO

	DO KY=1,NSTEP
	   YC(KY)=YMIN+REAL(KY)*STEPY
	   TSY(KY)=0.0
	   DO KX=1,NSTEP
	      TSY(KY)=TSY(KY)+TSO(KX,KY)
	   ENDDO
	ENDDO

! TINH CAC A(I)
	DO KX=1,NSTEP
	   A(KX)=0.0
	   DO KY=1,NSTEP
	      A(KX)=A(KX)+TSO(KX,KY)*YC(KY)
	   ENDDO
	ENDDO

! TINH D,B,C
	D=0.0
	B=0.0
	C=0.0
	DO KX=1,NSTEP
	   D=D+(A(KX)**2)/TSX(KX)
	ENDDO
	DO KY=1,NSTEP
	   B=B+TSY(KY)*YC(KY)
	   C=C+TSY(KY)*YC(KY)**2
	ENDDO

! TINH TY SO TUONG QUAN
	B=B**2/REAL(N)
	ETA=(D-B)/(C-B)
	TYSOTQ=ETA
	RETURN
	END

	SUBROUTINE REGRE(X,N,M,A,S,Q,U,RR,F)
! CHUONG TRINH NAY TINH CAC HE SO HOI QUI CUA PHUONG TRINH
! HOI QUI GIUA X1 VA CAC X2,...,XM VA XAC DINH CAC DAI LUONG
! CHO PHAN TICH PHUONG SAI
! DANG PHUONG TRINH HQ: X1=A1+A2*X2+...+AM*XM
!
! INPUT: + X MANG MOT CHIEU KICH THUOC N*M LUU TRU SO LIEU
!            QUAN TRAC CUA CAC X1..XM (LUU THEO COT)
!        + N DUNG LUONG MAU
!        + M SO BIEN (1 BIEN PHU THUOC VA M-1 BIEN DOC LAP)
!        + A MANG MOT CHIEU KICH THUOC M LUU CAC HE SO
!            A1, A2,..., AM
!        + S SAI SO CHUAN (CHUAN SAI THANG DU)
!        + Q TONG BINH PHUONG CAC BIEN SAI THANG DU
!        + U TONG BINH PHUONG CAC BIEN SAI HOI QUI
!        + RR HE SO TUONG QUAN BOI
!        + F BIEN CO PHAN BO FISHER
!
	DIMENSION X(N*M),A(M)
	DIMENSION R(M*M),RXX(M*M),RY(M),TB(M),RX(M*M)

	CALL COVAR(X,N,M,R,TB)
	L=0
	J=0
	DO I=M+1,M*M
	   IF (MOD(I,M).NE.1) THEN
	      L=L+1
	      RXX(L)=R(I)
	   ELSE
	      J=J+1
	      RY(J)=R(I)
	   ENDIF
	ENDDO
	RX=RXX
	CALL MINVERT(RX,M-1)
	CALL MULTIP(RX,RY,A,M-1,M-1,1)
	TMP=0.0
	DO J=1,M-1
	   TMP=TMP+A(J)*TB(J+1)
	ENDDO
	TMP=TB(1)-TMP
	DO J=M,2,-1
	   A(J)=A(J-1)
	ENDDO
	A(1)=TMP
	
	RX=R
	D=DET(RX,M)
	RX=R
	D11=PPDS(RX,M,1,1)
	TMP=D/D11

	RR=SQRT(1-TMP/R(1))

	Q=0.0
	DO I=1,N
	   YHQ=0.0
	   DO J=2,M
	      IJ=(J-1)*N+I
	      YHQ=YHQ+A(J)*X(IJ)
	   ENDDO
	   YHQ=YHQ+A(1)
	   TMP=X(I)-YHQ
	   Q=Q+TMP*TMP
	ENDDO
	U=REAL(N)*R(1)-Q
	F=U*REAL(N-M)/(Q*REAL(M-1))
	S=SQRT(Q/REAL(N-M))

	RETURN
	END

	SUBROUTINE MINVERT(AA,N)
! CHUONG TRINH NAY TINH MA TRAN NGHICH DAO CUA MA TRAN A
! INPUT: + AA MANG MOT CHIEU KICH THUOC N*N LUU THEO COT CUA A
!        + N SO HANG/COT CUA MA TRAN A
! OUTPUT: AA MA TRAN NGHICH DAO CUA A
!
	DIMENSION AA(N*N), A(N,N)

	K=0
	DO J=1,N
	   DO I=1,N
	      K=K+1
	      A(I,J)=AA(K)
	   ENDDO
	ENDDO
	M=N-1
	DO 10 I=1,N
	   R=A(I,1)
	   DO K=1,M
	      A(I,K)=A(I,K+1)/R
	   ENDDO
	   A(I,N)=1.0/R
	   DO 10 J=1,N
	      IF (I.NE.J) THEN
	         R=A(J,1)
	         DO K=1,M
	            A(J,K)=A(J,K+1)-R*A(I,K)
	         ENDDO
	         A(J,N)=-R*A(I,N)
	      ENDIF
  10  CONTINUE
	K=0
	DO J=1,N
	   DO I=1,N
	      K=K+1
	      AA(K)=A(I,J)
	   ENDDO
	ENDDO
	RETURN
	END

	SUBROUTINE MULTIP(A,B,C,N,M,L)
! CHUONG TRINH NAY TINH TICH CUA HAI MA TRAN
! INPUT: + A MANG MOT CHIEU KICH THUOC N*M (N HANG, M COT)
!            LUU MA TRAN A THEO COT
!        + B MANG MOT CHIEU KICH THUOC M*L (M HANG, L COT)
!            LUU MA TRAN B THEO COT
!        + N SO HANG CUA MA TRAN A
!        + M SO COT CUA A VA LA SO HANG CUA B
!        + L SO COT CUA B
! OUTPUT: + C MANG MOT CHIEU KICH THUOC N*L (N HANG, L COT)
!             LUU MA TRAN C THEO COT

	DIMENSION A(N*M),B(M*L),C(N*L)

	DO I=1,N
	DO K=1,L
	   IK=(K-1)*N+I
	   C(IK)=0.0
	   DO J=1,M
	      IJ=(J-1)*N+I   !(I-1)*M+J
	      JK=(K-1)*M+J
	      C(IK)=C(IK)+A(IJ)*B(JK)
	   ENDDO
	ENDDO
	ENDDO

	RETURN
	END
		
	FUNCTION PPDS(A,N,IR,JC)
!
! HAM NAY TINH PHAN PHU DAI SO CUA PHAN TU B(IR,JC)
! (HANG IR, COT JC) CUA MA TRAN B(NxN) MA CAC PHAN TU
! CUA NO DUOC LUU TRONG MANG A(N*N) THEO QUI CACH COT
! TRUOC DONG SAU
! INPUT: + MANG A(N*N) CHUA MA TRAN DAU VAO CUA B
!        + N KICH THUOC CUA B
!        + IR CHI SO HANG
!        + JC CHI SO COT
! OUTPUT: PHAN PHU DAI SO CUA PHAN TU HANG IR COT JC
!
	DIMENSION A(N*N), B(N,N)
	K=0
	DO J=1,N
	   DO I=1,N
	      K=K+1
	      B(I,J)=A(K)
	   ENDDO
	ENDDO
	
	K=0
	DO J=1,N
	   IF (J.NE.JC) THEN
	   DO I=1,N
	      IF (I.NE.IR) THEN
	         K=K+1
	         A(K)=B(I,J)
	      ENDIF
	   ENDDO
	   ENDIF
	ENDDO
	
	PPDS=DET(A,N-1)*(-1)**(IR+JC)
	RETURN
	END


	FUNCTION DET(X,N)
! HAM NAY TINH DINH THUC CUA MA TRAN A(N,N)
! INPUT: + MANG X DO DAI N*N CHUA CAC PHAN TU CUA MA TRAN A
!        + N KICH THUOC MA TRAN A
!   CHU Y: X(1)=A(1,1), X(2)=A(2,1),..., X(N)=A(N,1),X(N+1)=A(1,2),
!          X(N+2)=A(2,2),...
! OUTPUT: DINH THUC CUA A
!
	PARAMETER (EP=1.0E-6)
	DIMENSION X(N*N),A(N,N)
	K=0
	DO J=1,N
	   DO I=1,N
	      K=K+1
	      A(I,J)=X(K)
	   ENDDO
	ENDDO
	D=1.0
	N1=N-1
	DO 10 K=1,N1
	   AM=0.0
	   DO 11 I=K,N
	      T=A(I,K)
	      IF (ABS(T).GE.ABS(AM)) THEN
	         AM=T
	         J=I
	      ENDIF
  11     CONTINUE
         IF (ABS(AM).LE.EP) THEN
	      DT=0.0
	      DET=DT
	      RETURN
	   ELSE
	      IF (J.NE.K) THEN
	         D=-D
	         DO I=K,N
	            T=A(J,I)
	            A(J,I)=A(K,I)
	            A(K,I)=T
	         ENDDO
	      ENDIF
	   ENDIF
	   M=K+1
	   DO I=M,N
	      T=A(I,K)/AM
	      DO J=M,N
	         A(I,J)=A(I,J)-T*A(K,J)
		  ENDDO
	   ENDDO
	   D=D*A(K,K)
  10  CONTINUE
      DT=D*A(N,N)
	DET=DT
	RETURN
	END

	FUNCTION HSTQR(R,M,J,K)
! HAM NAY TINH HE SO TUONG QUAN RIENG GIUA BIEN XJ VA XK
! KHI XET DONG THOI M BIEN X1,X2,...,XM
! INPUT: + R MANG MOT CHIEU CHUA MA TRAN TUONG QUAN CHUAN HOA
!          (KICH THUOC M*M ==> R(M*M)
!        + M KICH THUOC MA TRAN TUONG QUAN, DONG THOI CUNG LA
!          SO BIEN DUOC XET
!        + J,K CHI SO BIEN DUOC TINH TUONG QUAN RIENG
! OUTPUT: HE SO TUONG QUAN RIENG GIUA BIEN THU J VA THU K
!
	DIMENSION R(M*M),RT(M*M)
	RT=R
	RJK=PPDS(RT,M,J,K)
	RT=R
	RJJ=PPDS(RT,M,J,J)
	RT=R
	RKK=PPDS(RT,M,K,K)
	HSTQR=-RJK/(SQRT(RJJ*RKK))
	RETURN
	END

	FUNCTION HSTQB(R,M,J)
! HAM NAY TINH HE SO TUONG QUAN BOI GIUA BIEN XJ VA TAP HOP
!        M-1 BIEN X1,X2,X(J-1),X(J+1),..,XM 
! INPUT: + R MANG MOT CHIEU CHUA MA TRAN TUONG QUAN CHUAN HOA
!          (KICH THUOC M*M ==> R(M*M)
!        + M KICH THUOC MA TRAN TUONG QUAN, DONG THOI CUNG LA
!          SO BIEN DUOC XET
!        + J CHI SO BIEN DUOC TINH TUONG QUAN BOI VOI TAP M-1 
!            BIEN CON LAI
! OUTPUT: HE SO TUONG QUAN BOI GIUA BIEN THU J VA CAC BIEN KHAC
!
	DIMENSION R(M*M),RT(M*M)
	RT=R
	RJJ=PPDS(RT,M,J,J)
	RT=R
	RR=DET(RT,M)
	HSTQB=SQRT(1-RR/RJJ)
	RETURN
	END