C @(#)mode.for 17.1.1.1 (ES0-DMD) 01/25/02 17:19:00 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.IDENTIFICATION C subroutine MODE version 3 820518 C A. Kruszewski ESO - Garching C.PURPOSE C calculates mode of data contained in array "A" C.ALGORITHM C mode is calculated with use of relation between mean, median C and mode, see: Kendall, M., Stuart, A. 1977, The Advanced Theory C of Statistics, 4-th edition, p. 40 C as this relation is valid only for mildly assymetric distributions C the data is repeatedly clipped around previous value of mode C first truncation is made around starting value of mean C if the dimension "N" of array "A" is greater than 5000 then only C every II-th element is used where II=N/5000+1 C C control array "CRMD" holds three parameters: C "CRMD(1)" introduces possibility of choice between mode and mean C subroutine returns mode when "CRMD(1)=3.0" C subroutine returns clipped mean when "CRMD(1)=0.0" C other values of "CRMD(1)" can be used when necessary C "CRMD(3)" is a clipping factor - data is clipped by "SIGMA*CRMD(3)" C from both sides of recent value of "AMODE" C iterations are terminated when the difference between two last C values of mode divided by last value of sigma is smaller than C "CRMD(2)" C maximum number of iterations is 10 C.INPUT/OUTPUT C input arguments: C A real array data array C N integer number of values in "A" C CRMD real array values of control parameters C output arguments: C AMODE real mode of array "A" C SIGMA real sigma of "A" after last clipping C----------------------------------------------------------------------- SUBROUTINE MODE(A,N,CRMD,AMODE,SIGM) C IMPLICIT NONE INTEGER N REAL A(N) REAL CRMD(3) REAL SIGM C INTEGER I, IC, II, IN, INN, L REAL ALC, AHC, AMEAN, AMED, AMODE REAL B(5000), C(5000), CFTR REAL DFTR, SCR, SIGMA, SFTR C C DFTR = CRMD(1) CFTR = CRMD(2) SFTR = CRMD(3) C C *** Define array B of dimension IN. IF (N.GT.5000) THEN II = (N-1) / 5000 + 1 IN = N / II INN = II * IN DO 10 I = II , INN , II B(I/II) = A(I) 10 CONTINUE ELSE DO 11 I = 1 , N B(I) = A(I) 11 CONTINUE IN = N ENDIF C C *** Calculate mean of non truncated data and C *** set it as a first approximation for mode. CALL MEAN(B, IN, AMODE, SIGMA) SIGM = SIGMA C C *** Set iteration control. IC = 9 C C *** Perform iterations. 12 CONTINUE C C *** Save recent value of mode. SCR = AMODE C C *** Set limits of truncation. ALC = AMODE - SFTR * SIGMA AHC = AMODE + SFTR * SIGMA C C *** Continue iterations. L=0 DO 13 I = 1 , IN IF ( B(I) .GT. ALC .AND. B(I) .LT. AHC ) THEN L = L + 1 C(L) = B(I) ENDIF 13 CONTINUE C C *** Stops iterations if number of data points C *** after recent truncation is less than 2. IF ( L .LT. 2 ) THEN IC = 0 ELSE CALL MEAN(C, L, AMEAN, SIGMA) C C *** Replaces mode with mean when DFTR.LE.0.1. IF ( DFTR .GT. 0.1 ) THEN CALL MEDIAN(C, L, ALC, AHC, AMED) IF ( L .GE. 10 ) THEN AMODE = DFTR * AMED - (DFTR-1.0) * AMEAN ELSE AMODE = AMED ENDIF ELSE AMODE=AMEAN ENDIF C C *** Compare two most recent values of mode and C *** stop iterations if they differ by less than C *** expected mean error of the mean of truncated C *** distribution multiplied by factor 'SGFACTOR'. SCR = MAX( ABS(SCR-AMODE) , ABS(AMODE-AMEAN) ) C C *** Stops iterations if sigma is not positive. IF (SIGMA.LE.0.0) THEN IC=0 ELSE SCR = SCR / ( SQRT( (SIGMA*SIGMA)/L ) ) IF ( SCR .LT. CFTR ) THEN IC = 0 ELSE IC = IC - 1 ENDIF ENDIF ENDIF IF (IC.GT.0) GOTO 12 C RETURN END