find_sub.f90

The module LSQ is for unconstrained linear least-squares fitting. It is based upon Applied Statisti

MODULE find_subsets

! A module of routines for finding and recording best-fitting subsets of
! regression variables

!   Version 1.10, 17 February 2004
!   Author: Alan Miller
!           formerly of CSIRO Division of Mathematical & Information Sciences
!   Phone: [phone]
!   e-mail: amiller @ bigpond.net.au
!   WWW-page: http://users.bigpond.net.au/amiller/

! 17 Feb 2004 Correction to subroutine EFROYM for the case in which all the
!             variables are selected. Thanks to David Jones.
! 12 Nov 1999 Made changes to routines exadd1 & add2 to prevent the calculation
!             of negative residual sums of squares which could occur in cases
!             in which the true RSS is zero.   Routine seq2 changed to avoid
!             cycling.
! 24 May 2000 Changed lsq_kind to dp (cosmetic change only)
! 4 June 2000 Added routine random_pick which picks a random set of variables.
! 29 Aug 2002 Set value of size in subroutine EFROYM when IER /= 0.

USE lsq

IMPLICIT NONE

INTEGER, SAVE                 :: max_size, nbest, lopt_dim1
REAL (dp), ALLOCATABLE, SAVE  :: bound(:), ress(:,:)
INTEGER, ALLOCATABLE, SAVE    :: lopt(:,:)


CONTAINS

SUBROUTINE init_subsets(nvar_max, fit_const, nvar)

INTEGER, INTENT(IN)           :: nvar_max
LOGICAL, INTENT(IN)           :: fit_const
INTEGER, INTENT(IN), OPTIONAL :: nvar

!     Local variables

INTEGER    :: i, ier
REAL (dp)  :: eps = 1.E-14
LOGICAL    :: lindep(ncol)

!     The LSQ module has probably already been initialized, but just in case ..

IF (.NOT. initialized) THEN
  IF (PRESENT(nvar)) CALL startup(nvar, fit_const)
END IF

IF (fit_const) THEN
  max_size = nvar_max + 1
ELSE
  max_size = nvar_max
END IF

lopt_dim1 = max_size * (max_size + 1) / 2
IF (ALLOCATED(bound)) DEALLOCATE(bound, ress, lopt)
ALLOCATE (bound(max_size), ress(max_size,nbest), lopt(lopt_dim1,nbest))

bound = HUGE(eps)
ress  = HUGE(eps)
lopt  = 0

CALL tolset(eps)
CALL sing(lindep, ier)

CALL ss()
DO i = 1, max_size
  CALL report(i, rss(i))
END DO

RETURN
END SUBROUTINE init_subsets


SUBROUTINE add1(first, last, ss, smax, jmax, ier)

! Calculate the reduction in residual sum of squares when one variable,
! selected from those in positions FIRST .. LAST, is added in position FIRST,
! given that the variables in positions 1 .. FIRST-1 (if any) are already
! included.

INTEGER, INTENT(IN)     :: first, last
INTEGER, INTENT(OUT)    :: jmax, ier
REAL (dp), INTENT(OUT)  :: ss(:), smax

!     Local variables

INTEGER    :: j, inc, pos, row, col
REAL (dp)  :: zero = 0.0_dp, diag, dy, ssqx, sxx(ncol), sxy(ncol)

!     Check call arguments

jmax = 0
smax = zero
ier = 0
IF (first > ncol) ier = 1
IF (last < first) ier = ier + 2
IF (first < 1) ier = ier + 4
IF (last > ncol) ier = ier + 8
IF (ier /= 0) RETURN

!     Accumulate sums of squares & products from row FIRST

sxx(first:last) = zero
sxy(first:last) = zero
inc = ncol - last
pos = row_ptr(first)
DO row = first, last
  diag = d(row)
  dy = diag * rhs(row)
  sxx(row) = sxx(row) + diag
  sxy(row) = sxy(row) + dy
  DO col = row+1, last
    sxx(col) = sxx(col) + diag * r(pos)**2
    sxy(col) = sxy(col) + dy * r(pos)
    pos = pos + 1
  END DO
  pos = pos + inc
END DO

!     Incremental sum of squares for a variable = SXY * SXY / SXX.
!     Calculate whenever sqrt(SXX) > TOL for that variable.

DO j = first, last
  ssqx = sxx(j)
  IF (SQRT(ssqx) > tol(j)) THEN
    ss(j) = sxy(j)**2 / sxx(j)
    IF (ss(j) > smax) THEN
      smax = ss(j)
      jmax = j
    END IF
  ELSE
    ss(j) = zero
  END IF
END DO

RETURN
END SUBROUTINE add1


SUBROUTINE add2(first, last, smax, j1, j2, ier)

!     Calculate the maximum reduction in residual sum of squares when 2
!     variables, selected from those in positions FIRST .. LAST, are
!     added, given that the variables in positions 1 .. FIRST-1 (if
!     any) are already included.    J1, J2 are the positions of the two
!     best variables.   N.B. J2 < J1.

INTEGER, INTENT(IN)     :: first, last
INTEGER, INTENT(OUT)    :: j1, j2, ier
REAL (dp), INTENT(OUT)  :: smax

!     Local variables

INTEGER    :: start, i1, i2, row, pos1, pos2, inc
REAL (dp)  :: zero = 0.0_dp, temp, det, two = 2.0, sxx(ncol), sxy(ncol), sx1x2

!     Check call arguments

smax = zero
j1 = 0
j2 = 0
ier = 0
IF (first > ncol) ier = 1
IF (last <= first) ier = ier + 2
IF (first < 1) ier = ier + 4
IF (last > ncol) ier = ier + 8
IF (ier /= 0) RETURN

start = row_ptr(first)

!     Cycle through all pairs of variables from those between FIRST & LAST.

DO i1 = first, last
  sxx(i1) = d(i1)
  sxy(i1) = d(i1) * rhs(i1)
  pos1 = start + i1 - first - 1
  DO row = first, i1-1
    temp = d(row) * r(pos1)
    sxx(i1) = sxx(i1) + temp*r(pos1)
    sxy(i1) = sxy(i1) + temp*rhs(row)
    pos1 = pos1 + ncol - row - 1
  END DO

  DO i2 = first, i1-1
    pos1 = start + i1 - first - 1
    pos2 = start + i2 - first - 1
    sx1x2 = zero
    DO row = first, i2-1
      sx1x2 = sx1x2 + d(row)*r(pos1)*r(pos2)
      inc = ncol - row - 1
      pos1 = pos1 + inc
      pos2 = pos2 + inc
    END DO
    sx1x2 = sx1x2 + d(i2)*r(pos1)

!     Calculate reduction in RSS for pair I1, I2.
!     The sum of squares & cross-products are in:
!              ( SXX(I1)  SX1X2   )      ( SXY(I1) )
!              ( SX1X2    SXX(I2) )      ( SXY(I2) )

    det = MAX( (sxx(i1) * sxx(i2) - sx1x2**2), zero)
    temp = SQRT(det)
    IF (temp < tol(i1)*SQRT(sxx(i2)) .OR.             &
        temp < tol(i2)*SQRT(sxx(i1))) CYCLE
    temp = ((sxx(i2)*sxy(i1) - two*sx1x2*sxy(i2))*sxy(i1) + sxx(i1)*sxy(i2)**2) &
           / det
    IF (temp > smax) THEN
      smax = temp
      j1 = i1
      j2 = i2
    END IF
  END DO ! i2 = first, i1-1
END DO   ! i1 = first, last

RETURN
END SUBROUTINE add2


SUBROUTINE bakwrd(first, last, ier)

!     Backward elimination from variables in positions FIRST .. LAST.
!     If FIRST > 1, variables in positions prior to this are forced in.
!     If LAST < ncol, variables in positions after this are forced out.
!     On exit, the array VORDER contains the numbers of the variables
!     in the order in which they were deleted.

INTEGER, INTENT(IN)  :: first, last
INTEGER, INTENT(OUT) :: ier

!     Local variables

INTEGER    :: pos, jmin, i
REAL (dp)  :: ss(last), smin

!     Check call arguments

ier = 0
IF (first >= ncol) ier = 1
IF (last <= 1) ier = ier + 2
IF (first < 1) ier = ier + 4
IF (last > ncol) ier = ier + 8
IF (ier /= 0) RETURN

!     For POS = LAST, ..., FIRST+1 call DROP1 to find best variable to
!     find which variable to drop next.

DO pos = last, first+1, -1
  CALL drop1(first, pos, ss, smin, jmin, ier)
  CALL exdrop1(first, pos, ss, smin, jmin)
  IF (jmin > 0 .AND. jmin < pos) THEN
    CALL vmove(jmin, pos, ier)
    IF (nbest > 0) THEN
      DO i = jmin, pos-1
        CALL report(i, rss(i))
      END DO
    END IF
  END IF
END DO

RETURN
END SUBROUTINE bakwrd


SUBROUTINE drop1(first, last, ss, smin, jmin, ier)

! Calculate the increase in the residual sum of squares when the variable in
! position J is dropped from the model (i.e. moved to position LAST),
! for J = FIRST, ..., LAST-1.

INTEGER, INTENT(IN)     :: first, last
INTEGER, INTENT(OUT)    :: jmin, ier
REAL (dp), INTENT(OUT)  :: ss(:), smin

!     Local variables

INTEGER    :: j, pos1, inc, pos, row, col, i
REAL (dp)  :: large = HUGE(1.0_dp), zero = 0.0_dp, d1, rhs1, d2, x, wk(last), &
              vsmall = TINY(1.0_dp)

!     Check call arguments

jmin = 0
smin = large
ier = 0
IF (first > ncol) ier = 1
IF (last < first) ier = ier + 2
IF (first < 1) ier = ier + 4
IF (last > ncol) ier = ier + 8
IF (ier /= 0) RETURN

!     POS1 = position of first element of row FIRST in r.

pos1 = row_ptr(first)
inc = ncol - last

!     Start of outer cycle for the variable to be dropped.

DO j = first, last
  d1 = d(j)
  IF (SQRT(d1) < tol(j)) THEN
    ss(j) = zero
    smin = zero
    jmin = j
    GO TO 50
  END IF
  rhs1 = rhs(j)
  IF (j == last) GO TO 40

!     Copy row J of R into WK.

  pos = pos1
  DO i = j+1, last
    wk(i) = r(pos)
    pos = pos + 1
  END DO
  pos = pos + inc

!     Lower the variable past each row.

  DO row = j+1, last
    x = wk(row)
    d2 = d(row)
    IF (ABS(x) * SQRT(d1) < tol(row) .OR. d2 < vsmall) THEN
      pos = pos + ncol - row
      CYCLE
    END IF
    d1 = d1 * d2 / (d2 + d1 * x**2)
    DO col = row+1, last
      wk(col) = wk(col) - x * r(pos)
      pos = pos + 1
    END DO
    rhs1 = rhs1 - x * rhs(row)
    pos = pos + inc
  END DO
  40 ss(j) = rhs1 * d1 * rhs1
  IF (ss(j) < smin) THEN
    jmin = j
    smin = ss(j)
  END IF

!     Update position of first element in row of r.

  50 IF (j < last) pos1 = pos1 + ncol - j
END DO

RETURN
END SUBROUTINE drop1


SUBROUTINE efroym(first, last, fin, fout, size, ier, lout)

!     Efroymson's stepwise regression from variables in positions FIRST,
!     ..., LAST.  If FIRST > 1, variables in positions prior to this are
!     forced in.  If LAST < ncol, variables in positions after this are
!     forced out.

!     A report is written to unit LOUT if LOUT >= 0.

INTEGER, INTENT(IN)    :: first, last, lout
INTEGER, INTENT(OUT)   :: size, ier
REAL (dp), INTENT(IN)  :: fin, fout

!     Local variables

INTEGER    :: jmax, jmin, i
REAL (dp)  :: one = 1.0, eps, zero = 0.0, ss(last), smax, base, var, f, smin

!     Check call arguments

ier = 0
IF (first >= ncol) ier = 1
IF (last <= 1) ier = ier + 2
IF (first < 1) ier = ier + 4
IF (last > ncol) ier = ier + 8
IF (fin < fout .OR. fin <= zero) ier = ier + 256
IF (nobs <= ncol) ier = ier + 512
IF (ier /= 0) THEN
  size = 0
  RETURN
END IF

!     EPS approximates the smallest quantity such that the calculated value of
!     (1 + EPS) is > 1.   It is used to test for a perfect fit (RSS = 0).

eps = EPSILON(one)

!     SIZE = number of variables in the current subset

size = first - 1

!     Find the best variable to add next

20 CALL add1(size+1, last, ss, smax, jmax, ier)
IF (nbest > 0) CALL exadd1(size+1, smax, jmax, ss, last)

!     Calculate 'F-to-enter' value

IF (size > 0) THEN
  base = rss(size)
ELSE
  base = rss(1) + ss(1)
END IF
var = (base - smax) / (nobs - size - 1)
IF (var < eps*base) THEN
  ier = -1
  f = zero
ELSE
  f = smax / var
END IF
IF (lout >= 0) WRITE(lout, 900) vorder(jmax), f
900 FORMAT(' Best variable to add:  ', i4, '  F-to-enter = ', f10.2)

!     Exit if F < FIN or IER < 0 (perfect fit)

IF (f < fin .OR. ier < 0) RETURN

!     Add the variable to the subset (in position FIRST).

IF (lout >= 0) WRITE(lout, '(50x, "Variable added")')
size = size + 1
IF (jmax > first) CALL vmove(jmax, first, ier)
DO i = first, MIN(jmax-1, max_size)
  CALL report(i, rss(i))
END DO

!     See whether a variable entered earlier can be deleted now.

30 IF (size <= first) GO TO 20
CALL drop1(first+1, size, ss, smin, jmin, ier)
CALL exdrop1(first+1, size, ss, smin, jmin)
var = rss(size) / (nobs - size)
f = smin / var
IF (lout >= 0) WRITE(lout, 910) vorder(jmin), f
910 FORMAT(' Best variable to drop: ', i4, '  F-to-drop  = ', f10.2)

IF (f < fout) THEN
  IF (lout >= 0) WRITE(lout, '(50x, "Variable dropped")')
  CALL vmove(jmin, size, ier)
  IF (nbest > 0) THEN
    DO i = jmin, size-1
      CALL report(i, rss(i))
    END DO
  END IF
  size = size - 1
  GO TO 30
END IF

IF (size >= last) RETURN
GO TO 20
END SUBROUTINE efroym


SUBROUTINE exadd1(ivar, smax, jmax, ss, last)

!     Update the NBEST subsets of IVAR variables found from a call
!     to subroutine ADD1.

INTEGER, INTENT(IN)    :: ivar, jmax, last
REAL (dp), INTENT(IN)  :: smax, ss(:)

!     Local variables

REAL (dp)  :: zero = 0.0_dp, ssbase, sm, temp, wk(last)
INTEGER    :: i, j, ltemp, jm

IF (jmax == 0) RETURN
IF (ivar <= 0) RETURN
IF (ivar > max_size) RETURN
ltemp = vorder(ivar)
jm = jmax
sm = smax
IF (ivar > 1) ssbase = rss(ivar-1)
IF (ivar == 1) ssbase = rss(ivar) + ss(1)
wk(ivar:last) = ss(ivar:last)

DO i = 1, nbest
  temp = MAX(ssbase - sm, zero)
  IF (temp >= bound(ivar)) EXIT
  vorder(ivar) = vorder(jm)
  IF (jm == ivar) vorder(ivar) = ltemp
  CALL report(ivar, temp)
  IF (i >= nbest) EXIT
  wk(jm) = zero
  sm = zero
  jm = 0
  DO j = ivar, last
    IF (wk(j) <= sm) CYCLE
    jm = j
    sm = wk(j)
  END DO
  IF (jm == 0) EXIT
END DO

!     Restore VORDER(IVAR)

vorder(ivar) = ltemp

RETURN
END SUBROUTINE exadd1


SUBROUTINE exdrop1(first, last, ss, smin, jmin)
! Record any new subsets of (LAST-1) variables found from a call to DROP1

INTEGER, INTENT(IN)    :: first, last, jmin
REAL (dp), INTENT(IN)  :: ss(:), smin

! Local variables
INTEGER    :: list(1:last), i