subset.f90

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

WRITE(11, '(a8, 2x, 4d13.4, f9.1)') yname, ymin, ymax, ymean, ystdev,   &
                                    (ymax - ymin)/ystdev
WRITE(11, *)

!     Test for singularities & form initial array (rss) of sums of squares
!     of residuals

ALLOCATE ( lindep(ncol) )
CALL sing(lindep, ier)
IF (ier .NE. 0) THEN
  WRITE(*, *)'** Rank deficiency in data = ', -ier, ' **'
  WRITE(11, *)'** Rank deficiency in data = ', -ier, ' **'
  DO i = i0+1, ncol
    IF (lindep(i)) THEN
      WRITE(*, '(1x, a, a8, a)') 'Variable', vname(vorder(i)-i0),  &
                   ' is linearly dependent on previous variables or is constant'
    END IF
  END DO
END IF
CALL ss()

!     Change extension to .red for output file.

ipos = INDEX(fname_dat, '.')
IF (ipos == 0) ipos = LEN_TRIM(fname_dat) + 1
fname_red = fname_dat(1:ipos) // 'red'

!     Find if a file with this name already exists

INQUIRE(FILE=fname_red, EXIST=ok)
IF (ok) THEN
  WRITE(*, *)'File: ', fname_red, ' already exists.  Overwrite? (Y/N): '
  READ(*, *) ans
  ok = (ans == 'Y' .OR. ans == 'y')
ELSE
  ok = .true.
END IF

!     Write R, EL, RHS & the residual sum of squares (sserr) to disk.

IF (ok) THEN
  OPEN(9, FILE=fname_red, ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
  WRITE(9) fname_dat, ypos, line1
  WRITE(9) nvar, ncol, nobs, r_dim, lsel
  IF (.NOT. fit_const) THEN
    WRITE(9) yname, vname(1:nvar)
  ELSE
    WRITE(9) yname, vname(0:nvar)
  END IF

  WRITE(9) r, d, rhs, rss, sserr
  CLOSE(9)
END IF

DEALLOCATE ( x, lindep, xmin, xmax, xmean, xstdev )

RETURN
END SUBROUTINE form_red_file



SUBROUTINE read_red_file()
!     Read R, EL, RHS & the residual sum of squares (sserr) from disk.

IMPLICIT NONE
INTEGER :: ncases

OPEN(9, FILE=fname_red, STATUS='OLD', ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
REWIND (9)
READ(9) fname_dat, ypos, line1
READ(9) nvar, ncol, ncases, r_dim, lsel
fit_const = (ncol > nvar)
CALL startup(nvar, fit_const)          ! From LSQ module
nobs = ncases

!     Allocate memory for array VNAME.

ALLOCATE ( vname(0:nvar) )

IF (.NOT. fit_const) THEN
  READ(9) yname, vname(1:nvar)
ELSE
  READ(9) yname, vname(0:nvar)
END IF

READ(9) r, d, rhs, rss, sserr
CLOSE (9)

WRITE(*, *)'File: ', fname_red, ' successfully read'

RETURN
END SUBROUTINE read_red_file



SUBROUTINE print_correlations()

!     Local variables

INTEGER              :: i1, i2, row, pos1, pos2
CHARACTER (LEN = 80) :: text

WRITE(*, *)'Do you want correlations amongst the predictor variables? (Y/N): '
READ(*, *) yesno
IF (yesno == 'y' .OR. yesno == 'Y') THEN
  IF (nv == 0) THEN
    WRITE(*, *)'Correlations amongst predictors'
    WRITE(11, *)'Correlations amongst predictors'
  ELSE
    WRITE(*, *)'Partial correlations amongst predictors'
    WRITE(*, '(1x, a/ (1x, 8a9))')'Variables forced in: ', vname(list)
    WRITE(11, *)'Partial correlations amongst predictors'
    WRITE(11, '(1x, a/ (1x, 8a9))')'Variables forced in: ', vname(list)
  END IF

!     i1 & i2 are the first & last column numbers
  i2 = in + 1
  DO
    i1 = i2 + 1
    i2 = MIN(i2 + 7, last)
    WRITE(*, '(11x, 7a9)') vname(vorder(i1:i2))
    WRITE(11, '(11x, 7a9)') vname(vorder(i1:i2))
    DO row = in+1, i2-1
      text = ' '
                     ! Find position in CORMAT of 1st off-diagonal correlation
                     ! in current row.
      pos = (row-in-1)*(last-in) - (row-in)*(row-in-1)/2 + 1
      pos2 = pos + i2 - row - 1
      WRITE(text, '(1x, a8)') vname(vorder(row))
      IF (row .GE. i1) THEN
        j = 20 + 9*(row-i1)
        pos1 = pos
      ELSE
        j = 11
        pos1 = pos + i1 - row - 1
      END IF
      WRITE(text(j:), '(7f9.4)') cormat(pos1:pos2)
      WRITE(*, '(a)') text
      WRITE(11, '(a)') text
    END DO
    IF (i2 .GE. last) EXIT
  END DO
END IF

IF (nv == 0) THEN
  WRITE(*, *)'Correlations with ', yname
  WRITE(11, *)'Correlations with ', yname
ELSE
  WRITE(*, *)'Partial correlations with ', yname
  WRITE(*, '(1x, a/ (1x, 8a9))')'Variables forced in: ', vname(list)
  WRITE(11, *)'Partial correlations with ', yname
  WRITE(11, '(1x, a/ (1x, 8a9))')'Variables forced in: ', vname(list)
END IF
i2 = in
DO
  i1 = i2 + 1
  i2 = MIN(i2 + 8, last)
  WRITE(*, '(2x, 8a9)') vname(vorder(i1:i2))
  WRITE(*, '(1x, 8f9.4)') ycorr(i1:i2)
  WRITE(11, '(2x, 8a9)') vname(vorder(i1:i2))
  WRITE(11, '(1x, 8f9.4)') ycorr(i1:i2)
  IF (i2 .GE. last) EXIT
END DO

RETURN
END SUBROUTINE print_correlations



SUBROUTINE current_status()
!     Print variable names & numbers and variables currently being forced
!     in or out (if any).

!     Local variables
CHARACTER (LEN=79) :: text
INTEGER            :: i

WRITE(*, 950) (i, vname(i), i=1, nvar)
950 FORMAT(6(' ', i3, ' ', a8))

text = 'Variables currently forced IN:  '
IF (first == i0+1) THEN
  text(33:36) = 'None'
ELSE
  DO i = i0+1, first-1
    WRITE(text(4*(i-i0)+29:4*(i-i0)+32), '(i4)') vorder(i)
  END DO
END IF
WRITE(*, '(" ", a)') text

text = 'Variables currently forced OUT: '
IF (last == ncol) THEN
  text(33:36) = 'None'
ELSE
  DO i = last+1, ncol
    WRITE(text(4*(i-last)+29:4*(i-last)+32), '(i4)') vorder(i)
  END DO
END IF
WRITE(*, '(1x, a)') text

RETURN
END SUBROUTINE current_status



SUBROUTINE f1max(ndf, nf, f1, f5, f10, ier)
! Calculates approximations to the 1%, 5% & 10% points of the distribution
! of the maximum of NF values of an F-ratio with 1 d.f. for the numerator
! and NDF d.f. for the denominator.
! An approximation is used to the values given in table 2 of:
!     Gilmour, S.G. (1996) `The interpretation of Mallows's Cp-statistic',
!     The Statistician, vol.45, pp.49-56

IMPLICIT NONE
INTEGER, INTENT(IN)  :: ndf, nf
REAL, INTENT(OUT)    :: f1, f5, f10
INTEGER, INTENT(OUT) :: ier

! Local variables
REAL :: a1(6) = (/ 1.67649, 6.94330,  1.22627, 0.25319,  0.06136, -2.41097 /)
REAL :: a5(6) = (/ 1.28152, 4.93268, -0.29583, 0.28518, -0.23566, -1.60581 /)
REAL :: a10(6) = (/1.06642, 3.96276, -0.62483, 0.30228, -0.52843, -1.04499 /)

REAL :: log_nf, one = 1.0

IF (ndf < 4 .OR. nf < 1) THEN
  ier = 1
  RETURN
END IF
ier = 0

log_nf = LOG(DBLE(nf))

f1 = one + EXP(a1(1) + (a1(3)/ndf + a1(2))/ndf + a1(4)*log_nf          &
               + (a1(6)/ndf + a1(5))/nf)

f5 = one + EXP(a5(1) + (a5(3)/ndf + a5(2))/ndf + a5(4)*log_nf          &
               + (a5(6)/ndf + a5(5))/nf)

f10 = one + EXP(a10(1) + (a10(3)/ndf + a10(2))/ndf + a10(4)*log_nf     &
                + (a10(6)/ndf + a10(5))/nf)

RETURN
END SUBROUTINE f1max



SUBROUTINE update(x, n, xmin, xmax, xmean, sxx)
! Update statistics for x
! n = observation number ( = 1 for the first call)
! (sxx is the sum of squares of deviations from the mean)

REAL, INTENT(IN)      :: x
INTEGER, INTENT(IN)   :: n
REAL, INTENT(IN OUT)  :: xmin, xmax, xmean, sxx

! Local variables
REAL  :: dev

IF (n == 1) THEN
  xmin = x
  xmax = x
  xmean = x
  sxx = 0.0
  RETURN
END IF

IF (x < xmin) THEN
  xmin = x
ELSE IF (x > xmax) THEN
  xmax = x
END IF
dev = x - xmean
xmean = xmean + dev/n
sxx = sxx + dev * (x - xmean)
RETURN
END SUBROUTINE update



SUBROUTINE calc_residuals(unit_data, unit_resid, fname_dat, nvar, vname,   &
                          ypos, nreq, beta, list, stdev, autoc1)

INTEGER, INTENT(IN)              :: unit_data, unit_resid, nvar, ypos, nreq,  &
                                    list(:)
CHARACTER (LEN = 40), INTENT(IN) :: fname_dat
CHARACTER (LEN = 8), INTENT(IN)  :: vname(0:)
REAL (dp), INTENT(IN)            :: beta(:), stdev
REAL (dp), INTENT(OUT)           :: autoc1

! Local variables
INTEGER    :: i, ier, iostatus
REAL (dp)  :: fit, lsq_resid, h, std_resid, last_resid = 0.0, x(0:nvar), &
              y, xrow(nreq), one = 1.0, zero = 0.0, sumsq

WRITE(unit_resid) 'Data file name: ', fname_dat
WRITE(unit_resid, 900) (list(i), vname(list(i)), i=1,nreq)
900 FORMAT('Variables in model:' / (5(i3, ' ', a8, ' | ')))
WRITE(unit_resid, 910) beta(1:nreq)
910 FORMAT('Regression coefficients used:' / (5g15.6))
WRITE(unit_resid) 'Actual Y  Fitted Y  LS-residual  Std-residual  Leverage'
x(0) = one
last_resid = zero
autoc1 = zero
sumsq = zero

DO
  READ(unit_data, *, IOSTAT=iostatus) x(1:ypos-1), y, x(ypos:nvar)
  IF (iostatus > 0) THEN
    CYCLE
  ELSE IF (iostatus < 0) THEN
    EXIT
  END IF
  fit = zero
  DO i = 1, nreq
    xrow(i) = x(list(i))
    fit = fit + beta(i) * xrow(i)
  END DO
  lsq_resid = y - fit
  CALL hdiag(xrow, nreq, h, ier)
  std_resid = lsq_resid / (stdev * SQRT(one - h))
  WRITE(unit_resid, 920) y, fit, lsq_resid, std_resid, h
  920 FORMAT(3g12.4, f9.2, f10.3)
  sumsq = sumsq + lsq_resid**2
  autoc1 = autoc1 + lsq_resid * last_resid
  last_resid = lsq_resid

END DO

autoc1 = autoc1 / sumsq
WRITE(*, '(a, f9.3)') 'Lag 1 auto-correlation of residuals = ', autoc1
WRITE(unit_resid, '(a, f9.3)') 'Lag 1 auto-correlation of residuals = ', autoc1

RETURN
END SUBROUTINE calc_residuals



SUBROUTINE cross_validation(unit_data, line1, ypos, fit_const, nvar, first,  &
			    last, search_method, criterion, nrepl, seed,   &
			    unit_rpt, msep, ier)
! Cross-validation routine excluding 10% at a time

! Search_method
! 1  Efroymson stepwise (F-to-enter = F-to-delete = 4.0)
! 2  Sequential replacement
! 3  Two-at-a-time sequential replacement
! 4  Best subsets (exhaustive search with branch-and-bound)

! Criterion (set to 5 for Efroymson search)
! 1  AIC
! 2  BIC
! 3  Mallows Cp
! 4  Hannan-Quinn
! 5  F-ratio = 4.0

! Other arguments:
! unit_data Unit number from which to read the data
! line1     Number of the first line of data in the input file (in case file
!           contains variable names & other header information)
! ypos      Position of the dependent variable in each line of data
! fit_const .true. if a constant is to be included in the model
! nvar      Number of predictor variables, excluding any constant
! first     Position of the first variable available for selection.   Any
!           variables in earlier positions are to be forced into all subsets.
! last      Position of the last variable available for selection.   Variables
!           to be forced out of subsets (if any) should be after position last.
! nrepl     Number of complete replications of 10 subsets of 10% of the data
! seed()    Array of random number seeds (dimension depends upon compiler)
! unit_rpt  Unit number for output of report
! msep      Average mean squared error of the predictions
! ier       Error indicator
!           = 0 if no error detected

! This version - 15 August 2002

INTEGER, INTENT(IN)      :: unit_data, line1, ypos, nvar, first, last,  &
                            search_method, nrepl, seed(:), unit_rpt
INTEGER, INTENT(IN OUT)  :: criterion
LOGICAL, INTENT(IN)      :: fit_const
INTEGER, INTENT(OUT)     :: ier
REAL (dp), INTENT(OUT)   :: msep

! Local variables
INTEGER              :: replicate, i, nobs_full, percentile, i1, i2, case,   &
                        iostatus, num_seeds, nvar_max, nsize, ipos, lout = 6, j
INTEGER, ALLOCATABLE :: order(:), seeds(:), list(:), vorder_cpy(:),  &
                        init_vorder(:)
REAL (dp)            :: sumsq_ls, total_LS, zero = 0.0_dp, minus1 = -1.0_dp, y,   &
			weight, fin = 4.0_dp, fout = 4.0_dp, variance, one = 1.0_dp, &
                        fit_LS, penalty, crit_val, min_crit_val
REAL (dp), ALLOCATABLE :: x(:), beta_LS(:)

IF (search_method == 1) criterion = 5

ALLOCATE( order(nobs), x(0:nvar), init_vorder(ncol) )
nobs_full = nobs
msep = zero
init_vorder = vorder

ALLOCATE( list(1:nvar) )
DO i = 1, nvar
  list(i) = i
END DO

CALL RANDOM_SEED(size=num_seeds)
ALLOCATE( seeds(num_seeds) )
seeds = seed(1:num_seeds)
CALL RANDOM_SEED(put=seeds)

! Return the QR factorization back to the order of variables in the data set.
DO i = 1, nvar
  CALL reordr(list, i, 2, ier)
END DO

DO replicate = 1, nrepl
! Choose a random permutation of the integers 1, 2, ..., nobs

  DO i = 1, nobs_full
    order(i) = i
  END DO

  CALL ranord(order, nobs_full)

  total_LS = zero

! ------------------- Cycle through subsets of the data ----------------

! Delete 10% of the observations using array `order'

  DO percentile = 1, 10
    IF (percentile == 1) THEN
      i1 = 1
      i2 = 0.1 * REAL(nobs_full) + 0.5
    ELSE
      i1 = i2 + 1
      IF (percentile == 10) THEN
	i2 = nobs_full
      ELSE
        i2 = 0.1 * REAL(nobs_full) * percentile + 0.5
      END IF
    END IF