subset.f90

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

            READ(10, *, IOSTAT=iostatus) y, x(1:nvar)
          ELSE
            READ(10, *, IOSTAT=iostatus) x(1:ypos-1), y, x(ypos:nvar)
          END IF

          IF (iostatus > 0) CYCLE                   ! Error in data
          IF (iostatus < 0) EXIT                    ! End of file

          e = y
          DO j = 1, i
            xcopy(j) = x(vorder(j))
            e = e - beta(j) * xcopy(j)
          END DO
          CALL hdiag(xcopy, i, h, ier)
          press = press + (e/(one - h))**2
        END DO
        WRITE(*, '(i4, 4x, 2g14.6, 2x, 10i4/ (38x, 10i4))')    &
                                   i-i0, press, press/nobs, vorder(1:i)
        WRITE(11, '(i4, 4x, 2g14.6, 2x, 10i4/ (38x, 10i4))')    &
                                   i-i0, press, press/nobs, vorder(1:i)
      END DO
      DEALLOCATE( x, beta, xcopy )

    CASE ('q', 'Q')          ! Quit this menu
      EXIT

    CASE DEFAULT
      WRITE(*, *) bel, '** Unknown option - ', ans, '  Try again! **'

  END SELECT

END DO

WRITE(*, *) 'Do you want to use cross-validation? (Y/N): '
READ(*, '(a)') ans
IF (ans == 'Y' .OR. ans == 'y') THEN
  DO
    WRITE(*, *) 'Choose search method'
    WRITE(*, *) '1. Stepwise (Efroymson)      2. Sequential replacement'
    WRITE(*, *) '3. 2-at-a-time replacement   4. Best subsets (exhaustive search)'
    WRITE(*, *) 'Enter number of your choice: '
    READ(*, *) search_method
    IF (search_method < 1 .OR. search_method > 4) THEN
      CYCLE
    ELSE
      EXIT
    END IF
  END DO

  IF (search_method > 1) THEN
    DO
      WRITE(*, *) 'Choose criterion for deciding size of subset'
      WRITE(*, *) '1. AIC (Akaike)           2. BIC (Bayesian)'
      WRITE(*, *) '3. Mallows Cp             4. Hannan-Quinn'
      WRITE(*, *) '5. F-ratio = 4.0'
      WRITE(*, *) 'Enter number of your choice: '
      READ(*, *) criterion
      IF (criterion < 1 .OR. criterion > 5) THEN
        CYCLE
      ELSE
        EXIT
      END IF
    END DO
  ELSE
    criterion = 5
  END IF

  WRITE(*, *) 'How many complete sets of 10 x 10% cross-validations: '
  READ(*, *) nrepl

  CALL RANDOM_SEED(size=i)
  WRITE(*, '(1x, a, i4, a)') 'Enter ', i, ' integers as random number seed(s): '
  ALLOCATE( seed(i) )
  READ(*, *) seed

  INQUIRE(10, OPENED=OK)
  IF (.NOT. OK) OPEN(10, FILE=fname_dat, STATUS='OLD')

  WRITE(11, '(/ a)') 'Using 10% cross-validation'
  WRITE(11, '("Search method ", a, "  Stopping criterion ", a/)')   &
                      method(search_method), crit_name(criterion)
  WRITE(11, *) 'Random number seeds: ', seed
  CALL cross_validation(10, line1, ypos, fit_const, nvar, first, last,      &
                        search_method, criterion, nrepl, seed, 11, msep, ier)
END IF   ! Want cross-validation?
STOP

CONTAINS


SUBROUTINE start()

!     This is the starting routine for the SUBSETS package of programs.
!     If a QR reduction has not been formed, then it forms the
!     upper-triangular Banachiewicz factorization of the input data.
!     Free-format input is assumed, i.e. with data fields separated by
!     spaces, CR's, tabs or commas.   N.B. Some Fortran compilers will
!     not accept tabs and/or commas as delimiters.

!     Latest revision - 10 May 2001

CHARACTER (LEN=20)    :: response
REAL (dp)             :: eps
REAL (dp), PARAMETER  :: pt0001 = 0.0001
INTEGER               :: i

!     Does the user already have a QR factorization?

WRITE(*, '(a)', ADVANCE='NO') ' Do you want to use an existing QR-factorization? (Y/N): '
READ(*, '(a)') ans
DO
  SELECT CASE (ans)
    CASE ('y', 'Y')                              ! Use existing .red file
      CALL get_file_name(fname_red, fname_red, 'red')
      OPEN(9, FILE=fname_red, STATUS='OLD', ACCESS='SEQUENTIAL',          &
           FORM='UNFORMATTED')
      CALL read_red_file()
      EXIT
    CASE ('n', 'N')
      CALL get_file_name(fname_dat, fname_dat, 'dat')
      CALL form_red_file()
      EXIT
    CASE DEFAULT
      WRITE(*, *) bel, '** You must enter Y or N **'
      CYCLE
  END SELECT
END DO

!     Set tolerances and test for singularities

eps = EPSILON(pt0001) ** 0.667
WRITE(*, '(a, g12.3)') ' Default tolerance for singularity test = ', eps
WRITE(*, *) 'This may be too small'
WRITE(*, *) 'e.g. if data recorded to say 5 significant decimals, use 1.E-5'
WRITE(*, '(a)', ADVANCE='NO')  &
            ' Enter tolerance or press RETURN to accept default: '
READ(*, '(a)') response
IF (LEN_TRIM(response) /= 0) THEN
  READ(response, '(f20.0)') eps
  eps = MAX( ABS(eps), 1.E-14)
  eps = MIN( eps, pt0001)
  WRITE(11, '(a, e12.3)') 'Tolerance used in singularity testing = ', eps
END IF
CALL tolset(eps)
IF (ALLOCATED(lindep)) DEALLOCATE(lindep)
IF (fit_const) THEN
  ALLOCATE ( lindep(0:nvar) )
ELSE
  ALLOCATE ( lindep(1:nvar) )
END IF
CALL sing(lindep, ier)
rank_deficit = -ier
IF (ANY(lindep)) THEN
  WRITE(*, '(i5, a)') -ier, ' singularities detected in predictor variables'
  WRITE(*, *) 'These variables are linearly related to earlier ones:'
  WRITE(11, '(i5, a)') -ier, ' singularities detected in predictor variables'
  WRITE(11, *) 'These variables are linearly related to earlier ones:'
  DO i = 1, nvar
    IF (lindep(i)) THEN
      WRITE(*, *) vname(i)
      WRITE(11, *) vname(i)
    END IF
  END DO
  WRITE(*, *)
  WRITE(11, *)
END IF

!     Write brief data summary to the screen and to the report file.

WRITE(*, *) 'QR-factorization is in file: ', fname_red
WRITE(11, *) 'QR-factorization is in file: ', fname_red

WRITE(*, '(" Dependent variable is: ", a8)') yname
WRITE(11, '(" Dependent variable is: ", a8)') yname
WRITE(*, 900) nvar, nobs, sserr
WRITE(11, 900) nvar, nobs, sserr
900 FORMAT(' No. of variables = ', i4, '     No. of observations = ', i5,   &
           / ' Residual sum of squares = ', g13.5)
IF (fit_const) THEN
  WRITE(*, *)'All models include a constant (intercept) term'
  WRITE(11, *)'All models include a constant (intercept) term'
  WRITE(*, 960) vname(0:nvar)
  WRITE(11, 960) vname(0:nvar)
  960 FORMAT(' Variable names: '/ (8(' ', a8)))
ELSE
  WRITE(*, *)'Note: Models do not include an intercept'
  WRITE(11, *)'Note: Models do not include an intercept'
  WRITE(*, 960) vname(1:nvar)
  WRITE(11, 960) vname(1:nvar)
END IF

IF (nobs > ncol) THEN
  ndf = nobs - ncol + rank_deficit
  WRITE(*, 910) sserr / ndf, ndf
  WRITE(11, 910) sserr / ndf, ndf
  910 FORMAT(' Residual variance estimate = ', g13.5, ' with ', i5,    &
             ' deg. of freedom'/)
END IF

!     Set default values for first & last
!     first = first variable which may an be omitted from subsets
!     last  = last variable which may be included in subsets

IF (fit_const) THEN
  first = 2
ELSE
  first = 1
END IF
last = ncol

RETURN
END SUBROUTINE start



SUBROUTINE get_file_name(fname1, fname2, extn)
!     Ask for details of the file name.   Add extension, if needed.

CHARACTER (LEN = 40), INTENT(OUT) :: fname1, fname2
CHARACTER (LEN =  3), INTENT(IN)  :: extn

!     Local variables

INTEGER              :: ipos, int_value(8)
LOGICAL              :: ok
CHARACTER (LEN = 3)  :: month(12) = (/'Jan', 'Feb', 'Mar', 'Apr', 'May',   &
                        'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'/)

DO
  WRITE(*, 900)
  900 FORMAT(' Name of data file = ? ')
  READ(*, *) fname1
  fname2 = fname1

!     Add extension if none has been entered,
!     detected by the lack of a '.'

  IF (INDEX(fname1, '.') .EQ. 0) THEN
    ipos = INDEX(fname1, ' ')
    fname2 = fname1(1:ipos-1) // '.' // extn
  END IF

!     Check that file exists.

  INQUIRE(FILE=fname2, EXIST=ok)
  IF (.NOT. ok) THEN
    WRITE(*, 910) bel, fname2
    910 FORMAT(' ', a, ' ** File not found - ', a, ' **')
    CYCLE
  ELSE
    EXIT
  END IF
END DO

!     Open report file and write time and date

ipos = INDEX(fname2, '.')
fname_rpt = fname2(1:ipos) // 'rpt'
INQUIRE(FILE=fname_rpt, EXIST=ok)
IF (.NOT. ok) THEN
  OPEN(11, FILE=fname_rpt)
ELSE
  DO
    WRITE(*, *) bel, 'File: ', fname_rpt, ' already exists'
    WRITE(*, *) 'Append (A) or overwrite (O)?: '
    READ(*, '(a)') ans
    SELECT CASE (ans)
      CASE ('a', 'A')        ! Append
        OPEN(11, FILE=fname_rpt, POSITION='APPEND')
        WRITE(11, '(//)')
        EXIT
      CASE ('o', 'O')        ! Overwrite
        OPEN(11, FILE=fname_rpt, STATUS='REPLACE')
        EXIT
      CASE DEFAULT
        WRITE(*, *) bel, 'Answer must be A or O -- try again!'
    END SELECT
  END DO
END IF

WRITE(11, '(a)') version
CALL DATE_AND_TIME(values = int_value)
WRITE(11, 950) int_value(5), int_value(6), int_value(7), int_value(3),    &
               month(int_value(2)), int_value(1)
950 FORMAT('Time ', i2, ':', i2.2, ':', i2.2, '     Date ', i2, ' ', a3,  &
           ' ', i4)
RETURN
END SUBROUTINE get_file_name



SUBROUTINE form_red_file()
!     Read in data from an ASCII file and form the QR-factorization

REAL (dp)               :: weight = 1.0_dp, y
REAL (dp), ALLOCATABLE  :: x(:)
CHARACTER (LEN = 79)    :: text
INTEGER                 :: ipos, i, iostatus, ier
LOGICAL                 :: ok
REAL, ALLOCATABLE       :: xmin(:), xmax(:), xmean(:), xstdev(:)
REAL                    :: ymin, ymax, ymean, ystdev

!     Display first part of file.

OPEN(10, FILE=fname_dat, STATUS='OLD')
WRITE(*, *)'Start of your data file follows'
DO i = 1, 12
  READ(10, '(a)') text
  WRITE(*, '(" ", a)') text
END DO
REWIND (10)

WRITE(*, 920, ADVANCE='NO')
920 FORMAT(' How many X-variables ? ')
READ(*, *) nvar
WRITE(*, 930, ADVANCE='NO')
930 FORMAT(' Do you want a constant in the model ? ')
READ(*, *) ans
fit_const = (ans .EQ. 'y' .OR. ans .EQ. 'Y')
IF (fit_const) THEN
  i0 = 1
ELSE
  i0 = 0
END IF

!     Allocate memory for arrays X and VNAME.

IF (fit_const) THEN
  ALLOCATE ( x(0:nvar), vname(0:nvar) )
ELSE
  ALLOCATE ( x(nvar), vname(nvar) )
END IF

!     Get position of dependant variable.

WRITE(*, '(a)', ADVANCE='NO') ' Is dependant variable at end ? (Y/N): '
READ(*, *) ans
IF (ans == 'Y' .OR. ans == 'y') THEN
  ypos = nvar + 1
ELSE
  WRITE(*, '(a)', ADVANCE='NO') ' Enter no. of position of dependant variable: '
  READ(*, *) ypos
  IF (ypos < 1) ypos = 1
  IF (ypos > nvar) ypos = nvar + 1
END IF

!     Enter variable names, read them from file, or set defaults.

IF (fit_const) vname(0) = 'Constant'
WRITE(*, *)'Are variable names in data file ? (Y/N): '
READ(*, *) ans
IF (ans == 'Y' .OR. ans == 'y') THEN
  WRITE(*, *)'Which line do names start on ? '
  READ(*, *) line1
  IF (line1 > 1) THEN
    DO i = 1, line1-1
      READ(10, *)
    END DO
  END IF
  IF (ypos > nvar) THEN
    READ(10, *) vname(1:nvar), yname
  ELSE IF (ypos == 1) THEN
    READ(10, *) yname, vname(1:nvar)
  ELSE
    READ(10, *) vname(1:ypos-1), yname, vname(ypos:nvar)
  END IF
  REWIND (10)
ELSE
  WRITE(*, *)'Do you want to name variables ? (Y/N): '
  READ(*, '(a)') ans
  IF (ans == 'Y' .OR. ans == 'y') THEN
    WRITE(*, *)'Variable names may contain up to 8 characters'
    WRITE(*, *)'Name for dependant (Y) variable = ? '
    READ(*, '(a)') yname
    DO i = 1, nvar
      WRITE(*, '(a)', ADVANCE='NO') ' Name for variable', i, ' = ? '
      READ(*, '(a)') vname(i)
    END DO
  ELSE
    DO i = 1, nvar
      IF (i < 10) THEN
        WRITE(vname(i), '("X(", i1, ")    ")') i
      ELSE IF (i < 100) THEN
        WRITE(vname(i), '("X(", i2, ")   ")') i
      ELSE IF (i < 1000) THEN
        WRITE(vname(i), '("X(", i3, ")  ")') i
      ELSE
        WRITE(vname(i), '("X(", i4, ") ")') i
      END IF
    END DO
    yname = 'Dept.var'
  END IF
END IF

WRITE(*, '(a)', ADVANCE='NO') ' Which line does the data start on ? '
READ(*, *) line1
IF (line1 > 1) THEN
  DO i = 1, line1-1
    READ(10, *)
  END DO
END IF

CALL startup(nvar, fit_const)          ! From LSQ module

!     Allocate arrays for storing descriptive statistics
ALLOCATE( xmin(nvar), xmax(nvar), xmean(nvar), xstdev(nvar) )

!     Read in data and form the upper-triangular factorization.

IF (fit_const) x(0) = one

!     Case is skipped if spurious characters are found (e.g. for missing values).

DO
  IF (ypos > nvar) THEN
    READ(10, *, IOSTAT=iostatus) x(1:nvar), y
  ELSE IF (ypos == 1) THEN
    READ(10, *, IOSTAT=iostatus) y, x(1:nvar)
  ELSE
    READ(10, *, IOSTAT=iostatus) x(1:ypos-1), y, x(ypos:nvar)
  END IF

  IF (iostatus > 0) CYCLE                   ! Error in data
  IF (iostatus < 0) EXIT                    ! End of file

! Update descriptive statistics
  DO i = 1, nvar
    CALL update(REAL(x(i)), nobs+1, xmin(i), xmax(i), xmean(i), xstdev(i))
  END DO
  CALL update(REAL(y), nobs+1, ymin, ymax, ymean, ystdev)

  CALL includ(weight, x, y)
END DO

!     Output descriptive statistics
WRITE(11, '("Variable      MinVal       MaxVal       Mean        Std.devn.   Range/stdev")')
xstdev = SQRT(xstdev/(nobs-1))
ystdev = SQRT(ystdev/(nobs-1))
DO i = 1, nvar
  WRITE(11, '(a8, 2x, 4d13.4, f9.1)') vname(i), xmin(i), xmax(i), xmean(i),  &
                                xstdev(i), (xmax(i) - xmin(i))/xstdev(i)
END DO