MODULE ALL_STATISTICS

CONTAINS

SUBROUTINE TEST1 (X,N,M)
REAL X(N,M),R(M,M)
X(:,:) = X(:,:)*2
print*,'CTC'
PRINT '(4F12.4)'	,X 

END SUBROUTINE



FUNCTION AVER(X,N)
!... HAM NAY TINH TRUNG BINH SO HOC CUA CHUOI
! INPUT: + X MANG DO DAI N CHUA SO LIEU QUAN TRAC
!        + N DUNG LUONG MAU
! OUTPUT: TRUNG BINH SO HOC CUA X
!
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 FUNCTION
!---------------------------------------------------------
FUNCTION STDEV (X,N)
!... HAM NAY TINH DO LECH CHUAN CUA CHUOI X
! INPUT: + X MANG DO DAI N CHUA SO LIEU QUAN TRAC
!        + N DUNG LUONG MAU
! OUTPUT: DO LECH CHUAN CUA X
!
INTEGER N, I
 REAL X(N)		! MANG CHUA SO LIEU BAN DAU
 REAL TMP, TMP1
 TMP = 0.0
 DO I = 1,N
    TMP = TMP + X(I)*X(I)
 END DO
 TMP1 = AVER(X,N)
 STDEV = SQRT(TMP/REAL(N) - TMP1*TMP1)
RETURN
END FUNCTION
!---------------------------------------------------------
SUBROUTINE XMAXMIN (X,N, AMAX, AMIN)
!... CT NAY TIM MAX, MIN CUA CHUOI X
! INPUT: + X MANG DO DAI N CHUA SO LIEU QUAN TRAC
!        + N DUNG LUONG MAU
! OUTPUT: MAX,MIN CUA X
!
INTEGER N,I,J
 REAL, INTENT(INOUT) :: X(N)	! MANG CHUA SO LIEU BAN DAU
 REAL AMAX, AMIN, TMP
 DO I=1,N-1
    DO J=I+1,N
	   IF (X(J) < X(I)) THEN
	      TMP = X(I)
		  X(I) = X(J)
		  X(J) = TMP
	   ENDIF
	ENDDO
 ENDDO
 AMAX = X(N)
 AMIN = X(1)
RETURN
END SUBROUTINE
!---------------------------------------------------------
FUNCTION P_LOWER(X,N,X0)
!... HAM NAY TINH UOC LUONG XAC SUAT CUA SU KIEN
!... P(X<X0)
! INPUT: + X MANG DO DAI N CHUA SO LIEU QUAN TRAC
!.       + N DUNG LUONG MAU
!        + X0 LA MOT GIA TRI NAO DO CUA X
! OUTPUT: XAC SUAT P(X<X0)
!
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
!---------------------------------------------------------
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
!---------------------------------------------------------

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
!---------------------------------------------------------

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
!---------------------------------------------------------

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	FUNCTION
!---------------------------------------------------------

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)
	DO
      I=I+1
      IF (I.GT.N) EXIT 
	  IF (X(I).NE.X(I-1)) THEN
	     N1=N1+1
		 X1(N1)=X(I)
	  ENDIF
	END DO
!  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	SUBROUTINE
!---------------------------------------------------------

FUNCTION PROB1(FR,XC,M,X0)
!
! HAM NAY TINH UOC LUONG XAC SUAT P(X<X0) TU HAM PHAN BO FR
! INPUT: + FR MANG CHUA TAN SUAT TICH LUY DO DAI M
!        + XC MANG DO DAI M CHUA TRI SO CUA X MA P(X<=XC(I))=FR(XC(I))
!        + M KICH THUOC CUA FR VA XC
!        + X0 GIA TRI UNG VOI XAC SUAT CAN TINH P(X<X0)
! OUTPUT: + XAC SUAT P(X<X0)
!
	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
! TIM NHOM CHUA X0
    J = 1
	DO
	  J = J + 1
      IF (XC(J).GE.X0) EXIT
	ENDDO
    PROB1=FR(J-1)+(X0-XC(J-1))*(FR(J)-FR(J-1))/(XC(J)-XC(J-1))
RETURN
END	FUNCTION
!---------------------------------------------------------

SUBROUTINE FREQ(X,N,XC,FR,F,NFR,M,XMIN,STEP,INTERVAL)
!
!  CHUONG TRINH NAY TINH TAN SO (NFR) VA TAN SUAT (FR)
!  CUA MANG X DO DAI N.
! INPUT:  + MANG X CHUA GIA TRI QUAN TRAC CUA BIEN NGAU NHIEN
!         + N LA DUNG LUONG MAU
!         + INTERVAL = 0 NEU LAY SO KHOANG CHIA NGAM DINH
!                    = 1 NEU TU CHON SO KHOANG CHIA 
! OUTPUT: + M LA SO KHOANG CHIA, KHI INTERVAL = 0 THI
!           M = 5*LOG10(N)
!         + XC MANG DO DAI M CHUA TRI SO TRUNG GIAN CUA CAC KHOANG
!         + NFR MANG DO DAI M CHUA TAN SO CAC NHOM
!         + FR MANG DO DAI M CHUA TAN SUAT TICH LUY CAC NHOM (HAM PHAN BO)
!         + F MANG DO DAI M CHUA MAT DO XAC SUAT CAC NHOM
!         + XMIN GIOI HAN DUOI CUA NHOM THU NHAT (<=MIN{X(I)}
!         + STEP BUOC CHIA
!
    REAL X(1),NFR(1),FR(1),XC(1),F(1)
	INTEGER N,INTERVAL,M
	REAL XMIN,XMAX,STEP
	CHARACTER ST*10
!	   
      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
!
!  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
!
      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
!
      STEP=(XMAX-XMIN)/REAL(M)
!  NEU MUON CHON BUOC TRON SO THI THEM DONG LENH SAU
!
	STEP=AINT(STEP)+1
!
! TINH XC
      XC(1)=XMIN+STEP/2.0
	DO J=2,M
	   XC(J)=XC(J-1)+STEP
	ENDDO
!  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	SUBROUTINE
!---------------------------------------------------------

FUNCTION PROB(FR,XCC,M,STEP,X0)
!
! HAM NAY TINH UOC LUONG XAC SUAT P(X<X0) TU HAM PHAN BO FR
! INPUT: + FR MANG CHUA TAN SUAT TICH LUY DO DAI M
!        + XC MANG CHUA TRI SO GIUA CAC NHOM DO DAI M
!        + M SO NHOM
!        + STEP CU LY NHOM (BUOC)
!        + X0 GIA TRI UNG VOI XAC SUAT CAN TINH P(X<X0)
! OUTPUT: + XAC SUAT P(X<X0)
!
	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
! 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	FUNCTION
!---------------------------------------------------------

FUNCTION X_ALFA(FR,XC,M,STEP,F0)
!
! HAM NAY TINH UOC LUONG TRI SO X0 KHI CHO TRUOC P(X<X0)=F0
! INPUT: + FR MANG CHUA TAN SUAT TICH LUY DO DAI M
!        + XC MANG CHUA TRI SO GIUA CAC NHOM DO DAI M
!        + M SO NHOM
!        + STEP CU LY NHOM (BUOC)
!        + F0=P(X<X0)=F(X0)
! OUTPUT: + TRI SO X0 UNG VOI P(X<X0)=F0
!
	DIMENSION FR(M),XC(M)
	REAL STEP, F0
	INTEGER M,J
!
	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	FUNCTION
!---------------------------------------------------------

SUBROUTINE FREQ_V(X,N,XC,FR,F,NFR,M,XMAX,XMIN,STEP)
!
!  CHUONG TRINH NAY TINH TAN SO (NFR) VA TAN SUAT (FR)
!  CUA MANG X DO DAI N.
! INPUT:  + MANG X CHUA GIA TRI QUAN TRAC CUA BIEN NGAU NHIEN
!         + N LA DUNG LUONG MAU
!         + XMAX,XMIN,STEP: MAX, MIN VA CU LY NHOM
!                    
! OUTPUT: + M LA SO KHOANG CHIA
!         + XC MANG DO DAI M CHUA TRI SO TRUNG GIAN CUA CAC KHOANG
!         + NFR MANG DO DAI M CHUA TAN SO CAC NHOM
!         + FR MANG DO DAI M CHUA TAN SUAT TICH LUY CAC NHOM 
!           (HAM PHAN BO)
!         + F MANG DO DAI M CHUA MAT DO XAC SUAT CAC NHOM
!
    REAL X(1),NFR(1),FR(1),XC(1),F(1)
	INTEGER N,M
	REAL XMIN,XMAX,STEP
	M=INT((XMAX-XMIN)/STEP)+1
!  TINH XC
      XC(1)=XMIN+STEP/2.0
	DO J=2,M
	   XC(J)=XC(J-1)+STEP
	ENDDO
!  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 SUBROUTINE
!---------------------------------------------------------

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,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	FUNCTION
!---------------------------------------------------------

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	SUBROUTINE
!---------------------------------------------------------

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
!---------------------------------------------------------


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
!---------------------------------------------------------

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
!---------------------------------------------------------

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	FUNCTION
!---------------------------------------------------------

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
!---------------------------------------------------------

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 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))
    ENDDO
RETURN
END	SUBROUTINE
!---------------------------------------------------------

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
!---------------------------------------------------------

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	SUBROUTINE
!---------------------------------------------------------

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,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
!---------------------------------------------------------

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
	DO
      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.LT.EP) EXIT
	  S=S1
	 ENDDO
     ANORMD=S1
RETURN
END	FUNCTION
!---------------------------------------------------------

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
    DO
 	 P0=ANORMD(0.0,1.0,A,B)
	 SS=(ABS(P0-AP)/AP)
	 IF (SS.LT.EPS) EXIT
     IF (AP.LT.P0) THEN
       C=B
       B=(A+B)/2.0
     ELSE IF (AP.GT.P0) THEN
       B=(C+B)/2.0
     ENDIF
	ENDDO
	A=(C+B)/2.0
	IF (IFLAG.EQ.1) THEN
 	   ANINV=-A*XICMA+MUY
	ELSE
       ANINV=A*XICMA+MUY
	ENDIF
RETURN
END	FUNCTION
!---------------------------------------------------------

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

	DO
 	P0=CHIDIST(B,N)
	SS=ABS((P0-AP)/AP)
	SS1=ABS((B-C)/B)
	IF (SS.LT.EPS.AND.SS1.LT.EPS) EXIT
	   IF (AP.GT.P0) THEN
	      C=B
	      B=(A+B)/2.0
	   ELSE IF (AP.LT.P0) THEN
	      B=(C+B)/2.0
	   ENDIF
	ENDDO
	A=(C+B)/2.0
    CHINV=A
RETURN
END	FUNCTION
!---------------------------------------------------------

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
!--------------------------------------------------------

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 FUNCTION
!---------------------------------------------------------

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 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
      enddo
      pause 'a too large, ITMAX too small in gcf'
1     gammcf=exp(-x+a*log(x)-gln)*h
return
END	SUBROUTINE
!---------------------------------------------------------

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 n=1,ITMAX
        ap=ap+1.
        del=del*x/ap
        sum=sum+del
        if(abs(del).lt.abs(sum)*EPS)goto 1
      enddo
      pause 'a too large, ITMAX too small in gser'
1     gamser=sum*exp(-x+a*log(x)-gln)
return
END	SUBROUTINE
!---------------------------------------------------------

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 j=1,6
        y=y+1.d0
        ser=ser+cof(j)/y
      enddo
      gammln=tmp+log(stp*ser/x)
return
END	FUNCTION
!---------------------------------------------------------

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

	DO
 	P0=TDIST(B,N)
	SS=ABS((P0-AP)/AP)
	SS1=ABS((B-C)/B)
	IF (SS.LT.EPS.AND.SS1.LT.EPS) EXIT
	   IF (AP.GT.P0) THEN
	      C=B
		  B=(A+B)/2.0
	   ELSE IF (AP.LT.P0) THEN
	      B=(C+B)/2.0
	   ENDIF
	ENDDO
	A=(C+B)/2.0
    TINV=A
RETURN
END	FUNCTION
!---------------------------------------------------------

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
!---------------------------------------------------------

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
!---------------------------------------------------------

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 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
      enddo
      pause 'a or b too big, or MAXIT too small in betacf'
1     betacf=h
return
END	FUNCTION
!---------------------------------------------------------

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

	DO
 	P0=FDIST(B,N1,N2)
	SS=ABS((P0-AP)/AP)
	SS1=ABS((B-C)/B)
	IF (SS.LT.EPS.AND.SS1.LT.EPS) EXIT
	   IF (AP.GT.P0) THEN
	      C=B
		  B=(A+B)/2.0
	   ELSE IF (AP.LT.P0) THEN
		  B=(C+B)/2.0
	   ENDIF
	ENDDO
	A=(C+B)/2.0
    FINV=A
RETURN
END	FUNCTION
!---------------------------------------------------------

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 FUNCTION
!---------------------------------------------------------
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(M*M) MANG MOT CHIEU LUU TRU MA TRAN
!            TUONG QUAN CHUAN HOA CUA M BIEN
!
REAL X(N*M),XX(N,M),R(M*M),TB(M),SX(M)
K=0
DO J=1,M
   TB(J)=AVER(X((J-1)*N+1:J*N),N)
   SX(J)=STDEV(X((J-1)*N+1:J*N),N)
   DO I=1,N
      K=K+1
      XX(I,J)=X(K)
   ENDDO
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
!*********************************************
SUBROUTINE CORRE_1(X,N,M,R)
! INPUT: + X(N,M) MANG HAI CHIEU LUU TRU SO LIEU
!        + N DUNG LUONG MAU
!        + M SO BIEN
! OUTPUT: + R(M,M) MANG HAI CHIEU LUU TRU MA TRAN
!            TUONG QUAN CHUAN HOA CUA M BIEN
!
REAL X(N,M), R(M,M),TB(M),SX(M)
DO J=1,M
   TB(J)=AVER (X(:,J),N)
   SX(J)=STDEV(X(:,J),N)
ENDDO
DO J=1,M
   DO K=J,M
      R(J,K)=HSTQ (X(:,J),X(:,K),N)
	  R(K,J)=R(J,K)
   ENDDO
ENDDO
RETURN
END SUBROUTINE
!---------------------------------------------------------

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
!---------------------------------------------------------

SUBROUTINE COVAR(X,N,M,R,TB)
! 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
!        + TB MANG MOT CHIEU KICH THUOC M CHUA TRUNG BINH CAC BIEN
!
	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 SUBROUTINE
!---------------------------------------------------------

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
   	   DO
       J=J+1
	   IF (X(I).EQ.XX(J)) THEN
	      U(I)=J
	      EXIT
	   ENDIF
	   IF (J.GE.N) EXIT
   	   ENDDO
   	ENDDO

	DO I=1,N
	   J=0
	   DO
       J=J+1
	   IF (Y(I).EQ.YY(J)) THEN
	      V(I)=J
	      EXIT
	   ENDIF
	   IF (J.GE.N) EXIT
       ENDDO
	ENDDO
	TMP=0.0
	DO I=1,N
	   TMP=TMP+(U(I)-V(I))**2
	ENDDO
	TQHANG=1.0-6.0*TMP/REAL(N*(N-1)*(N+1))
RETURN
END	FUNCTION
!---------------------------------------------------------

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	SUBROUTINE
!---------------------------------------------------------

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

FUNCTION HSTQ_1(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_1=RXY/SX/SY
RETURN
END FUNCTION
!---------------------------------------------------------

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	FUNCTION
!---------------------------------------------------------

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
!---------------------------------------------------------
SUBROUTINE GIAIHE_ND(A,B,X,N)
! CHUONG TRINH NAY GIAI HE PHUONG TRINH	DSTT Ax=b
! BANG PP NGHICH DAO MA TRAN
!
! INPUT: + A(N*N) MANG MOT CHIEU CHUA MA TRAN HE SO A
!            (LUU THEO COT)
!        + B(N) VECTO HE SO b
!        + N SO PHUONG TRINH
! OUTPUT: + X(N) VECTO NGHIEM
!
	REAL A(N*N),B(N),X(N)

	CALL MINVERT(A,N)
	CALL MULTIP(A,B,X,N,N,1)
RETURN
END	SUBROUTINE
!---------------------------------------------------------
SUBROUTINE GIAIHE_CRM(A,B,X,N)
! CHUONG TRINH NAY GIAI HE PHUONG TRINH	DSTT Ax=b
! BANG PP CRAMER
!
! INPUT: + A(N*N) MANG MOT CHIEU CHUA MA TRAN HE SO A
!            (LUU THEO COT)
!        + B(N) VECTO HE SO b
!        + N SO PHUONG TRINH
! OUTPUT: + X(N) VECTO NGHIEM
!
	REAL A(N*N),B(N),X(N), AA(N*N)
	REAL D, DI

	D = DET(A,N)
	DO I = 1,N
	   AA(:)=A(:)
	   AA((I-1)*N+1:I*N) = B(:)
	   DI = DET(AA,N)
	   X(I) = DI/D
	ENDDO
RETURN
END	SUBROUTINE
!-------------------------------------------------------
SUBROUTINE GIAIHE_GAUSE (AA,B,X,N)
integer N,i,j,k
REAL AA(N*N), B(N), X(N), A(N, N+1)

DO I=1,N
   A(:,I) = AA((I-1)*N+1:I*N)
ENDDO
A(:,N+1) = B(:)
DO I=1,N
	Tmp = A( I, I ) ! Sap xep lai de A (I,I) <> 0 
	IF (Tmp == 0) THEN      
	    K = I + 1   
	    DO WHILE (Tmp == 0 .AND. K < N)
	 	   Tmp = A( K, I )   
		   K = K + 1          
	    END DO 
      IF (Tmp == 0) THEN     
			PRINT*, "EXIT Here"
			STOP                      
		ELSE                          
		      X = A( I, 1:N+1 )             
		      K = K - 1          
		      A( I, 1:N+1 ) = A( K, 1:N+1 )          
		      A( K, 1:N+1 ) = X                
		END IF                         
	END IF
	Tmp = A(i,i)   ! Chuan hoa hang I 
  A(i,1:N+1) = A(i,1:N+1) / Tmp
  do i1=1,N    ! Khö cac phan tu ngoai d/cheo
     if (i1 /= i) then
       Tmp = A(i1,i)
  	  A(i1,1:N+1)=A(i1,1:N+1)-A(i,1:N+1)*Tmp
	   endif
  enddo
Enddo
X(1:N) = A(1:N,N+1)  ! Nghiem
END	SUBROUTINE

!---------------------------------------------------------

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 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 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
	   ENDDO
    ENDDO 
	K=0
	DO J=1,N
	   DO I=1,N
	      K=K+1
	      AA(K)=A(I,J)
	   ENDDO
	ENDDO
RETURN
END	SUBROUTINE
!--------------------------------------------------------
SUBROUTINE NDMT_GAUS (A,N)
REAL A(N,N), B(N,2*N), X(2*N)
integer i,j,k
REAL Tmp1, Tmp2

B = 0.   ! Khoi tao B
B( 1:N, 1:N ) = A(1:N,1:N) ! Nua dau cua B
DO I = 1, N     ! Duong cheo nua sau cua B
	 B( I, N+I ) = 1.
END DO
DO I = 1, N
	Tmp1 = B( I, I ) 
	IF (Tmp1 == 0) THEN ! Doi cho cac hang
	    K = I + 1   
	    DO WHILE (Tmp1 == 0 .AND. K <= N)
	 	    Tmp1 = B( K, I )   
		    K = K + 1          
	    END DO 
        IF (Tmp1 == 0) THEN     
			PRINT*, "EXIT Here"
			STOP                      
		ELSE                          
		      X = B( I, 1:2*N )             
		      K = K - 1          
		      B( I, 1:2*N ) = B( K, 1:2*N )          
		      B( K, 1:2*N ) = X                
		END IF                         
	END IF
! Chuan hoa hang I cho phan tu tren d/cheo
	B( I, 1:2*N ) = B( I, 1:2*N ) / Tmp1 
! Khu ngoai duong cheo nua dau cua B
   DO J = 1, N 
	   IF (J /= I) THEN
	      Tmp2 = B( J, I )
	      B( J, 1:2*N ) = B( J, 1:2*N ) &
              - B( I, 1:2*N ) * Tmp2
	   END IF
	 END DO
END DO
! Nghich dao cua A(1:N, 1:N)
A = B( 1:N, N+1:2*N )
! Nhan AI(1:N, 1:N) voi A(1:N,1:N) se nhan duoc nghiem
END	SUBROUTINE

!---------------------------------------------------------

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	SUBROUTINE
!---------------------------------------------------------
		
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
!---------------------------------------------------------

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 K=1,N1
	   AM=0.0
	   DO I=K,N
	      T=A(I,K)
	      IF (ABS(T).GE.ABS(AM)) THEN
	         AM=T
	         J=I
	      ENDIF
	   ENDDO
       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)
    ENDDO
    DT=D*A(N,N)
	DET=DT
RETURN
END	FUNCTION
!---------------------------------------------------------

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
!---------------------------------------------------------

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	FUNCTION
!---------------------------------------------------------

SUBROUTINE NGHIEM_CHIA_DOI (XL, XR, EPSILON, X0, ERROR, FF)
LOGICAL ERROR
REAL XL,XR,X0,EPSILON
REAL SS
EXTERNAL FF

ERROR = .FALSE.   ! CO NGHIEM TRONG KHOANG (XL,XR)
IF (FF(XL)==0.) THEN	! NGHIEM TAI DAU MUT TRAI
   X0 = XL
   RETURN
ELSE IF(FF(XR)==0.) THEN ! NGHIEM TAI DAU MUT PHAI
   X0 = XR
   RETURN
ELSE IF(XL==XR) THEN  ! HAI DAU MUT TRUNG NHAU VA
   ERROR = .TRUE.	  ! KHONG PHAI LA NGHIEM
   RETURN
ELSE IF(XL>XR) THEN	! DOI VI TRI HAI DAU MUT
   SS = XL
   XL = XR
   XR = SS
END IF
SS = ABS(XR-XL)
DO WHILE (SS >= EPSILON)
   X0 = (XR+XL)/2.
   FX0=FF(X0)
   FXL=FF(XL)
   FXR=FF(XR)
   IF (FX0 == 0) THEN	! NGHIEM DUNG TAI X0
      EXIT
   ELSE IF (FXL*FX0 < 0) THEN ! NGHIEM O NUA BEN TRAI
      XR = X0
   ELSE  IF (FXR*FX0 < 0) THEN ! NGHIEM O NUA BEN PHAI
      XL = X0
   END IF
   SS = ABS(XR-XL)  
ENDDO

END SUBROUTINE
!---------------------------------------------------------
FUNCTION SIMPSON (A, B, EP)
REAL A,B,DX,S1,S2,SS, EP
INTEGER N, I
!EXTERNAL FF

N = 0
S1 = 0.
DO
   N = N + 2
   DX = (B - A)/(2*N)
   S2 = 0.
   DO I = 1, N-1
	  S2 = S2 + FF(A + 2*I*DX)
   END DO
   SS = 0.
   DO I = 1, N
	  SS = SS + FF(A + (2*I-1)*DX)
   END DO
   S2 = DX/3*(FF(A) + FF(B) + 2*S2 + 4*SS)
   SS = ABS ((S2-S1)/S2)
   IF (SS < EP) EXIT
   S1 = S2
END DO
SIMPSON = S2
END	FUNCTION

!---------------------------------------------------------

FUNCTION FF(X)
REAL X
! PUT FUNCTION EXPRESSION HERE
FF = 1/SQRT(2*4*ATAN(1.))*EXP(-0.5*X*X)

END FUNCTION
!---------------------------------------------------------

SUBROUTINE DH1(U,UX,N,DX,SODO,CYC)
! CHUONG TRINH TINH DAO HAM BAC NHAT CUA HAM U(X)
! INPUT: + U(N) MANG 1 CHIEU CHUA GIA TRI CUA U(X)
!        + N: SO DIEM CO GIA TRI CUA U(X)
!        + DX: DO DAI KHOANG CHIA CUA X
!        + SODO: =1 SAI PHAN TIEN
!                =2 SAI PHAN LUI
!                =3 SAI PHAN TRUNG TAM DCX=2
!                =4 SAI PHAN TRUNG TAM DCX=4
!        + CYC: = .TRUE. NEU BIEN TUAN HOAN
!               = .FALSE. NEU BIEN KHONG TUAN HOAN
! OUTPUT: UX: DAO HAM CUA U(X)
!            SODO=1,CYC=.FALSE.: UX(1:N-1)
!            SODO=2,CYC=.FALSE.: UX(2:N)
!            SODO=3,CYC=.FALSE.: UX(2:N-1)
!            CYC=.TRUE.: UX(1:N)
INTEGER N,SODO, N1,N2,N3
LOGICAL CYC
REAL U(N),UX(N),C1,C2,DX2,DX4
SELECT CASE (SODO)
   CASE (1)
      DO I=1,N-1
		 UX(I) = (U(I+1)-U(I))/DX
	  ENDDO
	  IF (CYC) UX(N) = U(1)
   CASE (2)
      DO I=2,N
		 UX(I) = (U(I)-U(I-1))/DX
	  ENDDO
	  IF (CYC) UX(1) = U(N)
   CASE (3)
      DX2=2*DX
      DO I=2,N-1
		 UX(I) = (U(I+1)-U(I-1))/DX2
	  ENDDO
      UX(1) = (U(2)-U(1))/DX                                         
      UX(N) = (U(N)-U(N-1))/DX                                        
	  IF (CYC) THEN
	     UX(1) = (U(2)-U(N-1))/(2*DX)
		 UX(N) = UX(1)
	  END IF
   CASE (4)
      N1 = N-1                                                                    
      N2 = N-2                                                                    
      N3 = N-3                                                                    
      C1 = 4./3.                                                                  
      C2 = 1./3.                                                                  
      DX2 = 1./(2.*DX)                                                            
      DX4 = 1./(4.*DX)
      DO I = 3, N2                                                              
         UX(I) = C1*((U(I+1)-U(I-1))*DX2) &
               - C2*((U(I+2)-U(I-2))*DX4)                    
      ENDDO                                                                  
      IF (CYC) THEN
        UX(2)  = C1*((U(3)-U(1))*DX2)-C2*((U(4)-U(N-1))*DX4)                           
        UX(1)  = C1*((U(2)-U(N-1))*DX2)-C2*((U(3)-U(N-2))*DX4)                          
        UX(N)  = UX(1)                                                                 
        UX(N-1)= C1*((U(N)-U(N-2))*DX2)-C2*((U(2)-U(N-3))*DX4)
      END IF		                         
   END SELECT

END SUBROUTINE
!---------------------------------------------------------
SUBROUTINE DH2 (U,UX,N,DX,SODO)
! CHUONG TRINH TINH DAO HAM BAC HAI CUA HAM U(X)
! INPUT: + U(N) MANG 1 CHIEU CHUA GIA TRI CUA U(X)
!        + N: SO DIEM CO GIA TRI CUA U(X)
!        + DX: DO DAI KHOANG CHIA CUA X
!        + SODO: =1 DO CHINH XAC BAC 2
!                =2 DO CHINH XAC BAC 4
! OUTPUT: UX: DAO HAM CUA U(X)
!            SODO=1: UX(2:N-1)
!            SODO=2: UX(3:N-2)
INTEGER N,SODO
REAL U(N),UX(N), DX
REAL DX2,C15,C43,C12

DX2=DX*DX
SELECT CASE (SODO)
   CASE (1)
      DO I=2,N-1
	     UX(I) = (U(I+1)-U(I-1)-2.*U(I))/DX2
	  ENDDO
   CASE (2)
      C15=1./5.
	  C43=4./3.
	  C12=1./12.
	  DO I=3,N-2
		 UX(I)=(-C15*U(I)+C43*(U(I+1)-U(I-1)) &
		      -C12*(U(I+2)-U(I-2)))/DX2
	  ENDDO
   END SELECT
END SUBROUTINE
!--------------------------------------------

SUBROUTINE LAPLA52 (U,LAPU,H,N,M)
! CHUONG TRINH TINH LAPLAXIAN CUA U(X,Y)
! THEO SO DO 5 DIEM, DO CHINH XAC BAC 2
! INPUT: + U(N,M) GIA TRI CUA U
!        + N,M SO DIEM NUT THEO X VA Y
!        + H KHOANG CACH GIUA CAC DIEM NUT
! OUTPUT: LAPU(2:N-1,2:M-1) LAPLAXIAN CUA U
INTEGER N,M
REAL U(N,M),LAPU(N,M),H
INTEGER IP1, IM1, JP1, JM1

N1=N-1                                                                  
M1=M-1
DO I=2,N1
   DO J=2,M1
      IP1=I+1
	  IM1=I-1
	  JP1=J+1
	  JM1=J-1
      LAPU(I,J)=(U(IP1,J)+U(IM1,J)+U(I,JP1)+U(I,JM1) &
	            -4.*U(I,J))/H**2
   ENDDO
ENDDO
RETURN
END SUBROUTINE
!--------------------------------------------
SUBROUTINE LAPLA92 (U,LAPU,H,N,M)
! CHUONG TRINH TINH LAPLAXIAN CUA U(X,Y)
! THEO SO DO 9 DIEM, DO CHINH XAC BAC 2
! INPUT: + U(N,M) GIA TRI CUA U
!        + N,M SO DIEM NUT THEO X VA Y
!        + H KHOANG CACH GIUA CAC DIEM NUT
! OUTPUT: LAPU(2:N-1,2:M-1) LAPLAXIAN CUA U
INTEGER N,M
REAL U(N,M),LAPU(N,M),H
INTEGER IP1, IM1, JP1, JM1

N1=N-1                                                                  
M1=M-1
DO I=2,N1
   DO J=2,M1
      IP1=I+1
	  IM1=I-1
	  JP1=J+1
	  JM1=J-1
      LAPU(I,J)=(U(IP1,JP1)+U(IP1,JM1)+U(IM1,JP1) &
	            +U(IM1,JM1)+4.*(U(IP1,J)+U(IM1,J)+U(I,JP1) &
				+U(I,JM1))-20.*U(I,J))/(6.*H**2)
   ENDDO
ENDDO
RETURN
END SUBROUTINE



END MODULE ALL_STATISTICS