* This is a slightly modified version avoiding passing logicals from R SUBROUTINE wkuantile( FRACT, EL, X, W, IP, Q) *+ * Name: * wkuantile * Purpose: * Finds a quantile in a (possibly weighted) set of data. * Language: * Fortran 77 * Invocation: * CALL wkuantile( USEWT, INTERP, FRACT, EL, X, W, IP, Q) * Description: * The routine calculates the value of a specified quantile in a set * of data values, which may be weighted. In concept (although not * in practice) it sorts the supplied data values into ascending * order along with their associated positive weights, if supplied. * It then finds a quantile Q such that the sum of the weights * associated with all data values less than Q is a specified * fraction FRACT of the sum of all the weights supplied. If no * weights are supplied, then each data value is assigned unit * weight. There are two main applications of this algorithm: * * a) To find specified quantiles of a distribution of data values * for statistical purposes. In this case, the weights may * optionally be used to represent the number of times each data * value occurs. In such cases, it may be useful to regard the * distribution as continuous, and therefore to interpolate linearly * between data values when obtaining the result. * * b) Alternatively, the values may represent residuals from some * fitted function. In this case, by setting FRACT to 0.5, the * "weighted median residual" may be found. This has the property * that if it is subtracted from all the original residuals, then * the weighted sum of the absolute values of the corrected * residuals will be minimised. Thus, it may be used as the basis * for iteratively finding an `L1' fit. In such cases, the required * result will be equal to one of the data values (or may lie * mid-way between two of them) and interpolation between values is * not normally required. * Arguments: * USEWT = LOGICAL (Given) * Whether or not the data have associated weights. * INTERP = LOGICAL (Given) * Whether or not interpolation between data values should be * performed when obtaining the result. * FRACT = ? (Given) * The fraction specifying the required quantile, in the range * 0.0 to 1.0. * EL = INTEGER (Given) * Number of data values. * X( * ) = ? (Given) * Array of data values. * W( * ) = ? (Given) * Array of associated positive weights (if required). This * argument will only be referenced if USEWT is .TRUE.. * IP( EL ) = INTEGER (Given and Returned) * On entry, an array of pointers identifying which elements of X * (and W if supplied) are to be considered. On exit, these * pointers will have been permuted to access the specified data * elements in an order which is more nearly sorted than before * (although in general it will not represent a complete sort of * the data). * Q = ? (Returned) * The value of the requested quantile. * STATUS = INTEGER (Given and Returned) * The global status. * Notes: * - There are versions of this routine for processing both REAL * and DOUBLE PRECISION data; replace the "x" in the routine name by * R or D as appropriate. The types of the FACT, X, W and Q * arguments should match the routine being used. * - This routine is optimised for use when the number of data * values is large. In general, only a partial sort of the data will * be performed, so this routine will perform better than most other * methods of finding quantiles, which typically require a complete * sort. * - The order in which the input pointers are supplied in the * array IP is arbitrary, but there will often be an efficiency * advantage in supplying them so that they access the data in * nearly-sorted order. Thus, re-supplying the array of pointers * generated by a previous invocation of this routine (for the same * or similar data) may be worthwhile. * Algorithm: * This routine is a considerably enhanced Fortran implementation of * SELECT, which is algorithm no. 489 in the "Collected Algorithms * of the ACM", authors R.W. Floyd & R.L. Rivest (originally in * Algol). The enhancements include the use of a pointer array, the * ability to search for two ranked values simultaneously when * required (rather than the original single value) and the * iterative inversion of the cumulative weight versus rank function * in order to utilise weighted data. * * The resulting algorithm is rather complex in operation, and is * probably best approached by reading the original papers on the * SELECT algorithm (which is itself rather subtle, being both * recursive and iterative). The main additional complexities in * this implementation arise from: * * 1) The use of a recursive algorithm in a non-recursive language. * This is implemented in the usual way using a recursion stack, but * requires the ability to branch back from the end of a recursive * invocation of the algorithm into the middle of the calling * invocation. To avoid branching back into the range of a DO loop, * this requires that much of the looping be implemented using GO TO * statements. * * 2) Use of a pointer array to address the data being sorted. This * has the advantage that only one array need be sorted regardless * of whether or not weights are supplied. It also allows the data * and weight arrays to remain unmodified. Its implementatation is * best regarded is a simple substitution into the original * algorithm (which permuted the data themselves), although in * practice it leads to some additional difficulty in comprehension * and commenting. * * 3) Finding nearest-neighbours. The original SELECT algorithm was * intended to find only a single K'th smallest data value, but in * this implementation two values are generally required (lying * immediately above and below the required quantile). This is * handled by a modification to the method by which a new search * range is generated (depending on the results of the previous * iteration), together with a new termination condition. This * technique seeks to progressively label each partitioned element * as either "too small" or "too big", until there are no more * elements left. The last "too small" and "too big" elements are * then the ones required. * * 4) Separation of the `top level' invocation from other recursive * invocations. The top-level invocation of the algorithm behaves * differently to the recursive invocations of itself which it * subsequently makes. Only the top-level invocation is required to * find two array elements (see above), while the recursive * invocations merely have to return a single `good' partition value * for use in the next iteration of the next higher-level * invocation. They therefore continue to behave as in the original * SELECT algorithm. This difference is accommodated by varying the * method of selecting the next search range and the termination * condition according to the current depth of recursion. * * 5) Handling weights. If weights are supplied, then it is not * possible in advance to calculate the ranks of the data values * lying immediately above and below the required quantile. * Instead, these must be obtained by inverting the cumulative * weight versus rank function as iterations proceed. As a result, * a different top-level algorithm is needed in the case where * weights are supplied (recursive invocations made by this * top-level algorithm continue as normal). This is implemented by * having two code branches; one handles the `weight-less' and * recursive cases, while the other handles top-level invocations * for the weighted case alone. A return to the start of the * algorithm in order to perform another iteration can be made from * the end of either branch, with appropriate updating and * termination conditions being used in each case. * * 6) Inversion of the cumulative weight versus rank function. As * with the weight-less case, the top-level weighted algorithm seeks * to progressively label each partition value as either "too small" * or `too big' (until there are no more values to partition), this * assessment being made on the basis of the sum of weights lying on * each side. In addition, a new estimate of the target rank is * made at each iteration by linearly interpolating between the most * recent rank and sum-of-weights results. * Timing: * Details of the asymptotic time required to execute the original * SELECT algorithm are not altogether clear from the published * papers. It appears that this algorithm may have better average * performance than other methods and the time required may * approximate to EL * LOG( MIN( K, EL - K + 1 ) ) where K is the * rank of the largest data value which is smaller than the quantile * being sought. However, Sedgewick (see References) indicates that * such algorithms should, in general, complete in time proportional * to EL, so the above formula may be incorrect. When using * weighted data, the time will be multiplied by a further factor * reflecting the non-linearity of the cumulative weight versus rank * function and the difficulty of inverting it. * References: * - Comm. of the ACM, vol 18, no. 3 (March 1975), p165. * - Also see page 173. * - In addition, see the algorithm assessment by T. Brown in * Collected Algorithms of the ACM, (algorithm no. 489). * - Sedgwick, R., 1988, "Algorithms" (Addison-Wesley). * Copyright: * Copyright (C) 1992 Science & Engineering Research Council * Authors: * FR: R.W.Floyd and R.L Rivest (CACM) * RFWS: R.F. Warren-Smith (STARLINK, RAL) * DSB: David S. Berry (STARLINK) * RK: Roger Koenker(UIUC) * {enter_new_authors_here} * History: * ???-MAR-1975 (FR): * Original version in Algol. * 16-JUN-1989 (RFWS): * Translated to give a Fortran implementation. * 11-MAY-1992 (RFWS): * Considerably extended to handle both weighted and un-weighted * data. * 20-DEC-1994 (DSB): * Renamed from CCD1_QNTL to KPG1_QNTLx. * 29-JULY-2003 (RK): * adapted for use in R * 30-DEC-2020 (RK): * made USEWT and INTERP TRUE * {enter_further_changes_here} * Bugs: * {note_any_bugs_here} *- * Type definitions: IMPLICIT NONE ! No implicit typing * Arguments Given: LOGICAL USEWT PARAMETER (USEWT = .TRUE.) LOGICAL INTERP PARAMETER (INTERP = .TRUE.) DOUBLE PRECISION FRACT INTEGER EL * Arguments Given and Returned: DOUBLE PRECISION X( * ) DOUBLE PRECISION W( * ) INTEGER IP( EL ) * Arguments Returned: DOUBLE PRECISION Q * Local Constants: INTEGER MXSTK ! Size of stack for recursive calls * PARAMETER ( MXSTK = 100 ) PARAMETER ( MXSTK = 100 ) DOUBLE PRECISION HALF ! One half, in the required precision PARAMETER ( HALF = 1.0D0 / 2.0D0 ) * Local variables: DOUBLE PRECISION ALPHA ! Interpolation fraction DOUBLE PRECISION EXTRA ! Extra cumulative weight required DOUBLE PRECISION T ! Target value when partitioning DOUBLE PRECISION WL ! Sum of `left' weights for iteration DOUBLE PRECISION WLEFT ! Sum of `left' sorted weights DOUBLE PRECISION WMID ! Sum of weights of unsorted data DOUBLE PRECISION WR ! Sum of `right' weights for iteration DOUBLE PRECISION WRIGHT ! Sum of `right' sorted weights DOUBLE PRECISION WTARG ! Target weight for quantile INTEGER I ! General index for array elements INTEGER I1 ! Index of lower neighbouring element INTEGER I2 ! Index of upper neighbouring element INTEGER ITMP ! Temporary store for swapping pointers INTEGER J ! General index for array elements INTEGER K ! Target rank for SELECT algorithm INTEGER KZERO ! Initial estimate of target rank INTEGER L ! Left (first) array element to sort INTEGER LL ! Left (first) array element to sort INTEGER LSTK( MXSTK ) ! Stack of left array element indices INTEGER N ! Number of elements to be partitioned INTEGER R ! Right (last) array element to sort INTEGER RR ! Right (last) array element to sort INTEGER RSTK( MXSTK ) ! Stack of right array element indices INTEGER S ! Used for estimating sample size INTEGER SD ! Used for estimating sample size INTEGER ZD ! Used for estimating sample size INTEGER STK ! Stack pointer INTEGER G ! Used for estimating sample size DOUBLE PRECISION Z ! Used for estimating sample size *. * Initialise pointers to the left and right array elements to be * sorted. L = 1 R = EL * Initialise sums of weights for those elements which have already * been sorted (i.e. those lying to the left and right of the un-sorted * region). WLEFT = 0.0D0 WRIGHT = 0.0D0 * Calculate, or estimate (using linear interpolation), the rank of the * largest element such that the sum of weights for this and all * smaller elements will not exceed the required fraction of the total * sum of weights. Initialise the target rank to this value. KZERO = IDNINT( EL * MIN( MAX ( 0.0D0, FRACT ), 1.0D0 ) ) K = KZERO * Initialise the pointer to the recursive call stack. STK = 0 * A recursive invocation of the basic algorithm starts here. * ========================================================= 1 CONTINUE * Within each invocation, return to this point after each iteration, * until convergence is achieved. Confine the target rank for the next * iteration to the search range to be used. 2 CONTINUE ! Start of 'DO WHILE' loop K = MIN( MAX( L, K ), R ) * If there are enough elements still to be sorted to make it * worthwhile, then estimate the size of sample needed to get an * approximate indication of the target value (i.e. the data value with * the target rank). Do not perform this step if the stack will * overflow. (This is unlikely to happen, but if it does it simply * results in a slight loss of efficiency due to not reaching the * optimum depth of recursion. The results will still be OK.) IF ( ( ( R - L ) .GT. 60 ) .AND. ( STK .LT. MXSTK ) ) THEN N = R - L + 1 I = K - L + 1 * (The object here is to estimate the target value in such a way that * when IP(L:R) is permuted to partition X(IP(L:R)) about this value, * any error due to statistical uncertainty results in the true value * of the K'th ranked data value being pointed at by an element of IP * lying in the smaller part of the partitioned search range. The cost * of the next partitioning step is then minimised. The trick is to * balance the cost of estimating the target value - determined by the * sub-sample size - against the cost of getting it wrong. The original * paper should be consulted for details.) Z = DLOG( DBLE( N ) ) S = IDINT( HALF * DEXP( ( 2.0D0 / 3.0D0 ) * Z ) ) * This is a modification recommended by the algorithm assessor: * Original: ZD = INT( HALF * SQRT( Z * DBLE( S * ( N - S ) / N ) ) * : SIGN( 1, I - N / 2) ) * Modified to: G = S*(N-S)/N CALL INTPR("S",1,S,1) CALL INTPR("N",1,N,1) CALL INTPR("G",1,G,1) CALL DBLEPR("Z",1,Z,1) SD = INT( 0.1D0 * DSQRT( Z * REAL( S * ( N - S ) / N ) ) * : ( 2 * I / N - 1 ) ) CALL INTPR("SD",2,SD,1) * Set the bounds of the sub-sample of IP to be partitioned. LL = MAX( L, K - I * S / N + SD) RR = MIN( R, K + ( N - I ) * S / N + SD) * The basic algorithm now invokes itself recursively to get an * estimate of the target value from the smaller sub-sample it has * identified, before partitioning the entire sample it was given about * this value. Values which must be retained between recursive * invocations are pushed on to the stack. STK = STK + 1 LSTK( STK ) = L RSTK( STK ) = R * Branch back to the start of the algorithm with the appropriate * sub-sample defined. L = LL R = RR GO TO 1 * Return from a recursive invocation to this point. 3 CONTINUE END IF * Partition the data. * ================== * The code above implements the recursive invocation of the algorithm. * The following performs the actual sorting, based on an appropriate * iterative partitioning scheme. Two schemes are used; one for * unweighted data (and for recursive invocations, where the weights * are not relevant), and a second for the top-level invocation in the * case where weights are present (and the cumulative weight versus * rank function must be inverted at the same time). * Implement an "unweighted" sort. * ============================== * This is used for recursive invocations or when there are no weights * present; we are simply involved in sorting pointers to data values * without associated weights. Our aim is to get IP(K) to point at the * K'th ranked data value, and (in the case of a top-level invocation) * to get IP(K+1) to point at the (K+1)'th ranked value. IF ( ( .NOT. USEWT ) .OR. ( STK .GT. 0 ) ) THEN * We assume we have an estimate of the target value (i.e. that IP(K) * has been assigned a value so that it points to an approximation to * the K'th ranked data value). The IP array is then partitioned so * that all values smaller than X(IP(K)) are pointed at by elements of * IP lying to one side of IP(K) and all larger values are pointed at * by elements on the other side. We hope that the original value * X(IP(K)) will end up still being pointed at by IP(K) (in which case * the sort is complete), or by an element of IP pretty close to it * (preferably lying in the smaller part of the partitioned region). * The algorithm which follows is the basic QUICKSORT partitioning * sequence. Initialise the target value and the range of elements to * partition. T = X( IP( K ) ) I = L J = R * Initialise for partitioning, including a pointer to a `sentinel' * value at the end of the array. This allows the main part of the * partitioning algorithm to execute rapidly without having to check * against the array bounds. ITMP = IP( L ) IP( L ) = IP( K ) IP( K ) = ITMP IF ( X( IP( R ) ) .GT. T ) THEN ITMP = IP( R ) IP( R ) = IP( L ) IP( L ) = ITMP END IF * Move inwards from either end of the array region being partitioned * (I and J mark the array elements being considered at each end) until * there are no more elements to process. 4 CONTINUE ! Start of 'DO WHILE' loop IF ( I .LT. J ) THEN * Interchange pairs of elements which are out of order. ITMP = IP( I ) IP( I ) = IP( J ) IP( J ) = ITMP I = I + 1 J = J - 1 * Move the end pointers in towards the centre until an exchange of * values is indicated. 5 CONTINUE ! Start of 'DO WHILE' loop IF ( X( IP( I ) ) .LT. T ) THEN I = I + 1 GO TO 5 END IF 6 CONTINUE ! Start of 'DO WHILE' loop IF ( X( IP( J ) ) .GT. T ) THEN J = J - 1 GO TO 6 END IF GO TO 4 END IF * Tidy up, ensuring that IP(J) points at the value we have just * partitioned about. IF ( X( IP( L ) ) .EQ. T ) THEN ITMP = IP( L ) IP( L ) = IP( J ) IP( J ) = ITMP ELSE J = J + 1 ITMP = IP( J ) IP( J ) = IP( R ) IP( R ) = ITMP END IF * If this is a recursive invocation, then we are simply trying to * correctly set the value of IP(K). Update the left and right search * range so that this element lies between the closest elements of IP * so far correctly determined (of which the latest is IP(J)). Allow * the range to collapse to a single point if J equals K, in which case * the value of IP(K) has been found. IF ( STK .GT. 0 ) THEN IF ( J .LE. K ) L = J + 1 IF ( K .LE. J ) R = J - 1 * Continue iterating until only a single element (the K'th one) * remains to be found. Since its neighbours will have been found by * this point, it must now have the correct value. IF ( R .GT. L ) GO TO 2 * If this is a top-level invocation (not recursive), then we are * trying to correctly set two elements in the IP array (IP(KZERO) and * IP(KZERO+1)) so that they point at the correspondingly ranked data * values. Adjust the search range to exclude the element IP(J), which * has just been evaluated. Do this such that elements which point to * data values ranked KZERO and smaller get excluded to the left and * those pointing to values ranked (KZERO+1) and larger get excluded to * the right. (Note that we test against KZERO, because K may depart * from the required rank on the final iteration in order to stay * within the search range.) ELSE IF ( J .LE. KZERO ) THEN L = J + 1 ELSE R = J - 1 END IF * Iterations continue until no elements of IP remain in the search * range. The required elements (IP(KZERO) and IP(KZERO+1)) then * correspond with elements IP(L-1) and IP(R+1). IF ( R .GE. L ) GO TO 2 * When a solution has been found, obtain the indices of the data * elements lying immediately below and above the quantile and * calculate the fractional position of the target weight between these * elements. I1 = MAX( 1, L - 1 ) I2 = MIN( R + 1, EL ) ALPHA = ( EL * MIN( MAX( 0.0D0, FRACT ), 1.0D0 ) ) - : ( I1 - HALF ) ALPHA = MIN( MAX( 0.0D0, ALPHA ), 1.0D0 ) * To obtain the result, either interpolate between the adjacent * values... IF ( INTERP ) THEN Q = ( 1.0D0 - ALPHA ) * X( IP( I1 ) ) + : ALPHA * X( IP( I2 ) ) * ...or pick the nearer one, as required. ELSE IF ( ALPHA .LT. HALF ) THEN Q = X( IP( I1 ) ) ELSE IF ( ALPHA .EQ. HALF ) THEN Q = HALF * ( X( IP( I1 ) ) + X( IP ( I2 ) ) ) ELSE Q = X( IP( I2 ) ) END IF END IF END IF * Implement a top-level sort with weights. * ======================================== * If this is a top-level invocation of the algorithm (not recursive) * and weights have been supplied, then we must iteratively evaluate * two elements of the IP array (IP(K) and IP(K+1)), so that (a) they * point at the K'th and (K-1)'th ranked elements of X, and (b) these * values lie immediately below and above the required quantile. The * value of K is initially only a guess at its true final value, so * this must also be refined during the iterations. ELSE * We start by partitioning the IP array on the value X(IP(K)). As a * result of previous recursive invocations (if performed), this is the * current best estimate of the value of the K'th ranked data element * (using the current approximate value of K). The algorithm which * follows is the basic QUICKSORT partitioning sequence. Initialise the * target value, the range of elements to partition and sums for the * weights of elements which lie to the left and right of the * partition. T = X( IP( K ) ) I = L J = R WL = 0.0D0 WR = 0.0D0 * Initialise for partitioning, including a pointer to a `sentinel' * value at the end of the array. This allows the main part of the * partitioning algorithm to execute rapidly without having to check * aganst the array bounds. ITMP = IP( L ) IP( L ) = IP( K ) IP( K ) = ITMP IF ( X( IP( R ) ) .GT. T ) THEN ITMP = IP( R ) IP( R ) = IP( L ) IP( L ) = ITMP END IF * Move inwards from either end of the array region being partitioned * (I and J mark the array elements being considered at each end) until * there are no more elements to process. 7 CONTINUE ! Start of 'DO WHILE' loop IF ( I .LT. J ) THEN * Interchange pairs of elements which are out of order and sum their * weights. ITMP = IP( I ) IP( I ) = IP( J ) IP( J ) = ITMP WL = WL + W( IP( I ) ) WR = WR + W( IP( J ) ) I = I + 1 J = J - 1 * Move the end pointers in towards the centre until an exchange of * values is indicated, again summing the weights of elements on each * side of the partition as they are encountered. 8 CONTINUE ! Start of 'DO WHILE' loop IF ( X( IP( I ) ) .LT. T ) THEN WL = WL + W( IP( I ) ) I = I + 1 GO TO 8 END IF 9 CONTINUE ! Start of 'DO WHILE' loop IF ( X( IP( J ) ) .GT. T ) THEN WR = WR + W( IP( J ) ) J = J - 1 GO TO 9 END IF GO TO 7 END IF * Tidy up, ensuring that IP(J) points to the value we have just * partitioned about and that all weights have been correctly allocated * to left or right (the weight W(IP(J)) itself is allocated to the * left). IF ( I .EQ. J ) WL = WL + W( IP( I ) ) IF ( X( IP( L ) ) .EQ. T ) THEN ITMP = IP( L ) IP( L ) = IP( J ) IP( J ) = ITMP ELSE J = J + 1 ITMP = IP( J ) IP( J ) = IP( R ) IP( R ) = ITMP WL = WL + W( IP( J ) ) WR = WR - W( IP( J ) ) END IF * At this point, J contains an approximation to K (itself only an * approximation to its true final value, which is unknown) and * X(IP(J)) has the current partition value T (which is now known to be * the J'th ranked data value so that IP(J) has its correct value). * WLEFT contains the sum of weights for elements X(IP(1:L-1)), WL * contains the sum for elements X(IP(L:J)), WR contains the sum for * elements X(IP(J+1:R)) and WRIGHT contains the sum for elements * X(IP(R+1:EL)). Using these results, we must now form a new estimate * of K and update the search range for the next iteration. * Calculate the cumulative `target weight' which the sum of weights * for all elements smaller than the final result should have. WTARG = MIN( MAX( 0.0D0, FRACT ), 1.0D0 ) * : ( WLEFT + WL + WR + WRIGHT ) * Test whether the sum of elements lying to the `left' of IP(J), and * therefore smaller than X(IP(J)), plus the weight of element IP(J) * itself, lies above or below this target. In doing this, use only * half the weight of element IP(J). This reflects the fact that each * element is `centred'" in its weight bin. IF ( ( WLEFT + WL - HALF * W( IP( J ) ) ) .LT. WTARG ) THEN * If it lies below, then elements X(IP(1:J)) can now be regarded as * sorted (not literally, but in the sense that we need not consider * them again). Absorb the `left sum' of weights for this iteration * into the overall left sum of sorted weights and set the sum of * weights for the remaining unsorted elements (WMID) to the `right * sum' obtained on this iteration. Update the left boundary of the * unsorted region to exclude the newly-sorted elements. WLEFT = WLEFT + WL WMID = WR L = J + 1 * Otherwise, perform a similar operation to exclude elements to the * right, which can now be regarded as sorted. Ensure that the weight * of element IP(J) itself is transferred to the right and removed from * WMID during this process. ELSE WRIGHT = WRIGHT + WR + W( IP( J ) ) WMID = WL - W( IP( J ) ) R = J - 1 END IF * We now use linear interpolation to estimate (or guess) a better * value of K for the next iteration. (Note that before the next * iteration actually occurs, the new IP(K) will be updated via * recursive calls to the SELECT algorithm so that it points at an * improved estimate of the K'th ranked data value.) First calculate * the `extra cumulative weight' needed over and above the sum of * weights for all elements so far sorted to the left. EXTRA = WTARG - WLEFT * Adjust the weight values to remove half the weight of the largest * sorted left element and add half the weight of the smallest sorted * right element. This accounts for the elements being centred in their * weight bins. IF ( L .GT. 1 ) THEN EXTRA = EXTRA + HALF * W( IP( L - 1 ) ) WMID = WMID + HALF * W( IP( L - 1 ) ) END IF IF ( R .LT. EL ) THEN WMID = WMID + HALF * W( IP( R + 1 ) ) END IF * Make a new estimate of K by linear interpolation. Since the search * range has been reduced, this will be a more accurate estimate than * on the previous iteration. Allow an extra half element at each end * for accurate interpolation (weight bin centring again). Note that * the new value of K will be restricted to lie within the new search * range at the start of the next iteration. Use the central unsorted * element if there are no weights left in this region. IF ( WMID .NE. 0.0D0 ) THEN K = IDNINT( ( R - L + 2 ) * ( EXTRA / WMID ) ) + L - 1 ELSE K = ( L + R ) / 2 END IF * Iterations continue until there are no more elements of IP left to * sort. Since we repeatedly exclude elements pointing to values below * the quantile to the left, and exclude those pointing to values above * the quantile to the right, we eventually end up with the two * elements lying adjacent to the desired quantile being pointed at by * IP(L-1) and IP(R+1). IF ( R .GE. L ) GO TO 2 * When a solution has been found, obtain the indices of the data * elements lying immediately below and above the quantile and * calculate the fractional position of the target weight between the * cumulative sum of weights for each element. I1 = MAX( 1, L - 1 ) I2 = MIN( R + 1, EL ) ALPHA = ( WTARG - WLEFT + HALF * W( IP( I1 ) ) ) / : ( HALF * ( W( IP( I1 ) ) + W( IP( I2 ) ) ) ) ALPHA = MIN( MAX( 0.0D0, ALPHA ), 1.0D0 ) * To obtain the result, either interpolate between the adjacent * values... IF ( INTERP ) THEN Q = ( 1.0D0 - ALPHA ) * X( IP( I1 ) ) + : ALPHA * X( IP( I2 ) ) * ...or pick the nearer one, as required. ELSE IF ( ALPHA .LT. HALF ) THEN Q = X( IP( I1 ) ) ELSE IF ( ALPHA .EQ. HALF ) THEN Q = HALF * ( X( IP( I1 ) ) + X( IP ( I2 ) ) ) ELSE Q = X( IP( I2 ) ) END IF END IF END IF * Pop the stack and return from a recursive invocation of the * algorithm. This stage is omitted if the stack is empty, in which * case the return is from the top-level invocation so the sort is * complete. IF( STK .GT. 0 ) THEN L = LSTK( STK ) R = RSTK( STK ) STK = STK - 1 GO TO 3 END IF END