PROCEDURE FRSEL_FIT

1 Description

FRSEL_FIT is a Fortran 90 procedure for fitting a Fixed and Random effects SELectivity model, (Fryer, R.J. 1991. A model of between-haul variation in selectivity. ICES Journal of Marine Science, 48: 281-290).

2 Use

USE FRSEL_SETUP

USE FRSEL

CALL FRSEL_FIT(haul [, optional arguments])

Auxiliary routines required: NAG Fortran 77 Library, F01ABF, F03ABF.

3 Modules

FRSEL_FIT is contained in module FRSEL; the user need not set anything within this module. However, FRSEL_FIT requires module FRSEL_SETUP to set global parameters, and to define the data type that contains all haul information. FRSEL_SETUP must be of the following form:

MODULE FRSEL_SETUP ! specifies program parameters and sets-up haul information data type

IMPLICIT NONE

INTEGER, PARAMETER :: nalpha = 4 ! number of parameters to be estimated

INTEGER, PARAMETER :: nnu = 2 ! number of selectivity parameters for &

! each haul

INTEGER, PARAMETER :: wp = KIND(1.0d0) ! kind-value for double precision

TYPE HAUL_INFO ! derived type to contain haul information

REAL (wp), DIMENSION(nnu,nalpha) :: x ! design matrix

REAL (wp), DIMENSION(nnu) :: nu, b ! selectivity parameters and random effects

REAL (wp), DIMENSION(nnu,nnu) :: r, w ! within-haul and (inverse of) total variance

END TYPE HAUL_INFO

END MODULE FRSEL_SETUP

The user must specify:

nalpha	the number of parameters to be estimated,

nnu	the number of selectivity parameters estimated for each haul,

wp	the appropriate kind-value for the implementation of the NAG Fortran 77 library.

The type HAUL_INFO contains all haul information. It has five components:

%x	the nnu × nalpha design matrix X,

%nu	the nnu-vector of selectivity parameters ,

%b	the nnu-vector of random effects b,

%r	the nnu × nnu within-haul variance matrix R,

%w	the nnu × nnu matrix W.

4 Arguments

4.1 Mandatory argument

haul(:)- type(HAUL_INFO), intent(inout)

Input:	%x	the design matrix for each haul

%nu	the selectivity parameters for each haul.

Input/output:	%r	the within-haul variance matrix for each haul; only the lower triangle need be set; on exit, the upper triangle will be filled in.

Output:	%b	the random effect for each haul

%w	the W matrix for each haul.

The procedure will derive the number of hauls from the size of haul.

4.2 Optional arguments

Optional arguments must be supplied by keyword, not by position. They must all have exactly the dimensions described below.

mltype - character*1, intent(in)

Input: the type of likelihood estimation, 'm' for maximum likelihood or 'r' for residual maximum likelihood.

Constraint: must be lower case.

Default: mltype = 'r'

niter - integer, intent(in)

Input: the number of iterations.

Information: niter ≤ 1 gives the default number of iterations.

Default: niter = 1000.

tol - real (wp), intent(in)

Input: the convergence criterion; convergence is deemed to have occurred if the difference in (residual) log-likelihood between two successive iterations is less than or equal to tol.

Constraint: tol > 0.0_wp.

Default: tol = 10.0_wp**(-12).

monit - integer, intent(in)

Input: controls convergence monitoring by printing every monit iterations.

Information: monit ≤ 0 gives no monitoring information.

Default: no printing.

d_specified - logical, intent(in)

Input: indicates if the user has given a starting estimate of the between-haul variance D:

if d_specified = .true. the user must provide an initial estimate of D,

if d_specified = .false. the default value of d will be used.

Constraint: if d_specified = .true., d must be set.

Default: d_specified = .false..

d(nnu,nnu) - real (wp), intent(inout)

Input: the lower triangle of an initial estimate of the between-haul variance D.

Output: the estimate of D.

Constraint: d must be positive definite.

Default: the identity matrix.

log_lik - real (wp), intent(out)

Output: the maximised (residual) log-likelihood.

alpha(nalpha) - real (wp), intent(out)

Output: the (residual) maximum likelihood parameter estimates .

se(nalpha) - real (wp), intent(out)

Output: the standard errors of the (residual) maximum likelihood parameter estimates .

vcov(nalpha, nalpha) - real (wp), intent(out)

Output: the variance matrix of the (residual) maximum likelihood parameter estimates .

error_message - integer, intent(out)

Output: error_message = 0 unless the procedure detects an error.

Information: if error_message is not set, the procedure stops execution of the program if it detects an error; thus, a program will only terminates normally if no errors are detected.

5 Error indicators and warnings

error_message = 1

convergence not obtained in niter iterations; either increase niter, or start again with current value of d.

error_message = 2

mltype is not equal to either 'm' or 'r'.

error_message = 3

d is not positive-definite.

error_message = 4

d is not specified, even though d_specified = .true..

11 ≤ error_message ≤ 14

error calculating det( XTX) in NAG routine F03ABF, corresponding to NAG ifail = error_message - 10.

21 ≤ error_message ≤ 22

error calculating W = V-1 in NAG routine F01ABF corresponding to NAG ifail = error_message - 20.

31 ≤ error_message ≤ 34

error calculating det(V) in NAG routine F03ABF corresponding to NAG ifail = error_message - 30.

41 ≤ error_message ≤ 42

error calculating Var( ) in NAG routine F01ABF corresponding to NAG ifail = error_message - 40.

51 ≤ error_message ≤ 54

error calculating det( XTWX) in NAG routine F03ABF, corresponding to NAG ifail = error_message - 50.

6 Example

Taken from Fryer (1991). Fits the four parameter model described on p287 to the data in Table 2 by residual maximum likelihood. Note that the program results are not identical to those on p288 due to rounding errors in the printing of the covariance matrices in Table 2.

6.1 Program text

MODULE FRSEL_SETUP ! specifies program parameters and sets-up haul information data type

IMPLICIT NONE

INTEGER, PARAMETER :: nalpha = 4 ! number of parameters to be estimated

INTEGER, PARAMETER :: nnu = 2 ! number of selectivity parameters for &

! each haul

INTEGER, PARAMETER :: wp = KIND(1.0d0) ! kind-value for double precision

TYPE HAUL_INFO ! derived type to contain haul information

REAL (wp), DIMENSION(nnu,nalpha) :: x ! design matrix

REAL (wp), DIMENSION(nnu) :: nu, b ! selectivity parameters and random effects

REAL (wp), DIMENSION(nnu,nnu) :: r, w ! within-haul and (inverse of) total variance

END TYPE HAUL_INFO

END MODULE FRSEL_SETUP

INCLUDE 'FRSEL.F90'

PROGRAM EXAMPLE

USE FRSEL_SETUP

USE FRSEL

! example of the use of FRSEL1 taken from Fryer, 1991 (ICES Journal of Marine Science, &

! 48: 281-290

! the data are given in Table 2 (p 286) and the four parameter model described on p287 &

! is fitted: note that the parameter estimates differ from those on p288 because of &

! rounding errors in the printing of the covariance matrices in Table 2

IMPLICIT NONE

INTEGER, PARAMETER :: nhaul = 10

INTEGER :: i, j, k

REAL (wp) :: d(nnu,nnu), alpha(nalpha), se(nalpha)

TYPE (HAUL_INFO) :: haul(nhaul)

! read in parameter estimates and covariance matrices for each haul

OPEN (FILE = 'nu.dat', UNIT = 9, STATUS = 'old')

READ (9,*) ((haul(i)%nu(j), j = 1,nnu), ((haul(i)%r(j,k), k = 1,j), j = 1,nnu), i = 1,nhaul)

! construct x matrices

DO i = 1, nhaul

haul(i)%x = 0.0_wp

haul(i)%x(1,1) = 1.0_wp

haul(i)%x(2,2) = 1.0_wp

END DO

DO i = 7, nhaul

haul(i)%x(1,3) = 1.0_wp

haul(i)%x(2,4) = 1.0_wp

END DO

! fit fixed and random effects selectivity model

CALL FRSEL_FIT(haul, alpha = alpha, se = se, d = d)

! write output

WRITE(*,100) (alpha(j), se(j), j = 1, nalpha), ((d(j,k), k = 1, 2), j = 1, 2)

STOP

100 FORMAT(/,' alpha se',/,4(/,2f9.3),//,' d',/,2(/,2f9.4),//)

END PROGRAM EXAMPLE

6.2 Program data

-7.61 0.437 1.128 -0.0564 0.00285

-9.63 0.581 1.805 -0.0909 0.00461

-9.75 0.499 1.104 -0.0543 0.00271

-4.47 0.286 0.371 -0.0179 0.00089

-8.51 0.591 7.174 -0.3870 0.02106

-4.32 0.325 0.744 -0.0408 0.00221

-7.43 0.297 0.800 -0.0279 0.00099

-11.86 0.384 0.743 -0.0247 0.00084

-10.76 0.329 0.666 -0.0215 0.00071

-8.49 0.290 1.683 -0.0642 0.00254

6.3 Program results

alpha se

-6.672 0.954

0.415 0.039

-3.005 1.478

-0.088 0.057

d

4.1598 -0.1320

-0.1320 0.0057

MODULE FRSEL ! fits fixed and random effects selectivity model of Fryer (1991)

! programmer: Rob Fryer

! date: 21 October 1996

! language: Fortran 90 (standard)

! auxiliary routines: NAG F01ABF, F03ABF

! documentation: available from author

USE FRSEL_SETUP

CONTAINS

SUBROUTINE FRSEL_FIT(haul, mltype, niter, tol, monit, d_specified, d, log_lik, alpha, se, vcov, error_message)

IMPLICIT NONE

! Arguments

TYPE (HAUL_INFO), DIMENSION(:), INTENT(INOUT) :: haul ! haul information

CHARACTER*1, OPTIONAL, INTENT(IN) :: mltype ! type of likelihood: 'm' = maximum likelihood, or the default,

! 'r' = residual maximum likelihood

INTEGER, OPTIONAL, INTENT(IN) :: niter ! number of iterations: default 1000

REAL (wp), OPTIONAL, INTENT(IN) :: tol ! convergence tolerance: default 10**(-12)

INTEGER, OPTIONAL, INTENT(IN) :: monit ! controls monitoring, printing every monit iterations: default 0

LOGICAL, OPTIONAL, INTENT(IN) :: d_specified ! determines whether d is specified: default false

REAL (wp), OPTIONAL, INTENT(INOUT) :: d(nnu, nnu) ! between-haul variance: default identity matrix

REAL (wp), OPTIONAL, INTENT(OUT) :: log_lik ! (residual) log-likelihood

REAL (wp), OPTIONAL, INTENT(OUT) :: alpha(nalpha) ! parameters to be estimated

REAL (wp), OPTIONAL, INTENT(OUT) :: se(nalpha) ! standard errors

REAL (wp), OPTIONAL, INTENT(OUT) :: vcov(nalpha,nalpha) ! variance matrix

INTEGER, OPTIONAL, INTENT(OUT) :: error_message

! Local variables

INTEGER :: i, niter_wk

REAL (wp) :: tol_wk, d_wk(nnu,nnu), log_lik_wk, alpha_wk(nalpha), vcov_wk(nalpha,nalpha)

CHARACTER*1 :: mltype_wk

IF (PRESENT(error_message)) error_message = 0

CALL INITIALISE_VARIABLES(haul, mltype, niter, tol, d_specified, d, mltype_wk, niter_wk, tol_wk, d_wk, error_message)

IF (PRESENT(error_message) .and. error_message /= 0) RETURN

CALL EM_ALGORITHM(haul, mltype_wk, niter_wk, tol_wk, monit, d_wk, log_lik_wk, alpha_wk, vcov_wk, error_message)

IF (PRESENT(error_message) .and. error_message /= 0 .and. error_message /= 1) RETURN

IF (PRESENT(log_lik)) log_lik = log_lik_wk ! assign optional output variables given convergence, or if run out of

IF (PRESENT(alpha)) alpha = alpha_wk ! iterations

IF (PRESENT(se)) THEN

DO i = 1, nalpha

se(i) = sqrt(vcov_wk(i,i))

END DO

END IF

IF (PRESENT(d)) d = d_wk

IF (PRESENT(vcov)) vcov = vcov_wk

RETURN

CONTAINS

! ****************************************************************************************************************************

SUBROUTINE INITIALISE_VARIABLES(haul, mltype, niter, tol, d_specified, d, mltype_wk, niter_wk, tol_wk, d_wk, error_message)

IMPLICIT NONE

! Arguments

TYPE (HAUL_INFO), INTENT(INOUT), DIMENSION(:) :: haul

CHARACTER*1, OPTIONAL, INTENT(IN) :: mltype

INTEGER, OPTIONAL, INTENT(IN) :: niter

REAL (wp), OPTIONAL, INTENT(IN) :: tol

LOGICAL, OPTIONAL, INTENT(IN) :: d_specified

REAL (wp), OPTIONAL, INTENT(IN) :: d(nnu,nnu)

CHARACTER*1, INTENT(OUT) :: mltype_wk

REAL (wp), INTENT(OUT) :: tol_wk

INTEGER, INTENT(OUT) :: niter_wk

REAL (wp), INTENT(OUT) :: d_wk(nnu,nnu)

INTEGER, OPTIONAL, INTENT(INOUT) :: error_message

! Local variables

INTEGER :: i, j, ifail

REAL (wp) :: wk(nnu), det

DO i = 1, SIZE(haul) ! fill in haul variance matrices

DO j = 1, nnu

haul(i)%r(j,j:nnu) = haul(i)%r(j:nnu,j)

END DO

END DO

IF (PRESENT(mltype)) THEN ! initialise monitoring type

mltype_wk = mltype

ELSE

mltype_wk = 'r'

END IF

IF (mltype_wk .ne. 'm' .and. mltype_wk .ne. 'r') THEN

CALL ERROR_CHECK(2, error_message)

RETURN

END IF

IF (PRESENT(niter) .and. niter > 1) THEN ! initialise number of iterations

niter_wk = niter

ELSE

niter_wk = 1000

END IF

IF (PRESENT(tol)) THEN ! initialise convergence tolerance

tol_wk = tol

ELSE

tol_wk = 10.0_wp**(-12)

END IF

IF (PRESENT(d_specified) .and. d_specified) THEN ! initialise between-haul variance

IF (PRESENT(d)) THEN

d_wk = d

DO i = 1, nnu ! fill in upper triangle

d_wk(i,i:nnu) = d_wk(i:nnu,i)

END DO

ifail = 1

CALL F03ABF(d_wk, nnu, nnu, det, wk, ifail) ! check d_wk is positive-definite

IF (ifail /= 0) THEN

CALL ERROR_CHECK(3, error_message)

RETURN

END IF

ELSE

CALL ERROR_CHECK(4, error_message)

RETURN

END IF

ELSE

d_wk = 0.0_wp

DO i = 1, nnu

d_wk(i,i) = 1.0_wp

END DO

END IF

RETURN

END SUBROUTINE INITIALISE_VARIABLES

! ****************************************************************************************************************************

SUBROUTINE EM_ALGORITHM(haul, mltype, niter, tol, monit, d, log_lik, alpha, vcov, error_message)

IMPLICIT NONE

! Arguments

TYPE (HAUL_INFO), DIMENSION(:), INTENT(INOUT) :: haul

CHARACTER*1, INTENT(IN) :: mltype

INTEGER, INTENT(IN) :: niter

REAL (wp), INTENT(IN) :: tol

INTEGER, OPTIONAL, INTENT(IN) :: monit

REAL (wp), INTENT(INOUT) :: d(nnu,nnu)

REAL (wp), INTENT(OUT) :: log_lik

REAL (wp), INTENT(OUT) :: alpha(nalpha)

REAL (wp), INTENT(OUT) :: vcov(nalpha,nalpha)

INTEGER, OPTIONAL, INTENT(INOUT) :: error_message

! Local variables

INTEGER :: iiter

REAL (wp) :: log_lik_old, log_lik_1

CALL LOG_LIK_1_CALC(haul, mltype, log_lik_1, error_message) ! calculate constant term in log-likelihood

IF (PRESENT(error_message) .and. error_message /= 0) RETURN

DO iiter = 1, niter ! EM algorithm

IF (iiter > 1) log_lik_old = log_lik ! store old log-likelihood

log_lik = log_lik_1 ! initialise log-likelihood with constant

CALL W_CALC(haul, d, log_lik, error_message) ! calculate variance matrices w and increment log-likelihood

IF (PRESENT(error_message) .and. error_message /= 0) RETURN

CALL VCOV_CALC(haul, mltype, vcov, log_lik, error_message) ! calculate var(alpha) and increment likelihood

IF (PRESENT(error_message) .and. error_message /= 0) RETURN

CALL B_CALC(haul, vcov, d, alpha, log_lik) ! calculate alpha, b, and increment log-likelihood

CALL D_CALC(haul, vcov, mltype, d) ! calculate new value of d

log_lik = - log_lik / 2.0_wp ! convert to log-likelihood

IF (iiter > 1) THEN ! check for convergence

IF (log_lik - log_lik_old .le. tol) RETURN

END IF

IF (PRESENT(monit) .and. monit > 0 .and. MOD(iiter, monit) == 0) CALL MPRINT(log_lik, alpha, d)

END DO

CALL ERROR_CHECK(1, error_message)

RETURN

END SUBROUTINE EM_ALGORITHM

! ****************************************************************************************************************************

SUBROUTINE LOG_LIK_1_CALC(haul, mltype, log_lik_1, error_message)

! calculate constant term in (residual) log-likelihood

IMPLICIT NONE

! Arguments

TYPE (HAUL_INFO), DIMENSION(:), INTENT(IN) :: haul

CHARACTER*1, INTENT(IN) :: mltype

REAL (wp), INTENT(OUT) :: log_lik_1

INTEGER, OPTIONAL, INTENT(INOUT) :: error_message

! Local variables

INTEGER :: i, ifail, nhaul

REAL (wp) :: wk(nalpha, nalpha), det, wkspce(nalpha)

REAL (wp) :: pi = 3.14159265

nhaul = SIZE(haul)

IF (mltype == 'r') THEN

wk = 0.0_wp ! calculate sum(transpose(x) *+ x)

DO i = 1, nhaul

wk = wk + MATMUL(TRANSPOSE(haul(i)%x), haul(i)%x)

END DO

ifail = 1

CALL F03ABF(wk, nalpha, nalpha, det, wkspce, ifail) ! calculate det(sum)

IF (ifail .ne. 0) THEN

CALL ERROR_CHECK(ifail + 10, error_message)

RETURN

END IF

log_lik_1 = REAL(nnu * nhaul - nalpha, wp) * log(2.0_wp * pi) - log(det)

ELSE

log_lik_1 = REAL(nnu * nhaul, wp) * log(2.0_wp * pi)

END IF

RETURN

END SUBROUTINE LOG_LIK_1_CALC

! ****************************************************************************************************************************

SUBROUTINE W_CALC(haul, d, log_lik, error_message)

! calculate w matrices and increment log-likelihood

IMPLICIT NONE

! Arguments

TYPE (HAUL_INFO), DIMENSION(:), INTENT(INOUT) :: haul

REAL (wp), INTENT(IN) :: d(nnu,nnu)

REAL (wp), INTENT(INOUT) :: log_lik

INTEGER, OPTIONAL, INTENT(INOUT) :: error_message

! Local variables

INTEGER :: i, j, ifail

REAL (wp) :: v(nnu+1,nnu), wk(nnu), det

DO i = 1, SIZE(haul)

v(1:nnu,:) = haul(i)%r + d ! v = r + d

ifail = 1

CALL F01ABF(v, nnu+1, nnu, haul(i)%w, nnu, wk, ifail) ! w = inverse(v)

IF (ifail /= 0) THEN

CALL ERROR_CHECK(20+ifail, error_message)

RETURN

END IF

DO j = 1, nnu

haul(i)%w(j,j:nnu) = haul(i)%w(j:nnu,j) ! fill in upper triangle

END DO

ifail = 1 ! det(v)

CALL F03ABF(v(1:nnu,:), nnu, nnu, det, wk, ifail)

IF (ifail /= 0) THEN

CALL ERROR_CHECK(30+ifail, error_message)

RETURN

END IF

log_lik = log_lik + log(det) ! increment log-likelihood by sum(log(det(v)))

END DO

END SUBROUTINE W_CALC

! ***************************************************************************************************************************

SUBROUTINE VCOV_CALC(haul, mltype, vcov, log_lik, error_message)

! calculate variance(alpha) and increment log-likelihood further

IMPLICIT NONE

! Arguments

TYPE (HAUL_INFO), DIMENSION(:), INTENT(IN) :: haul

CHARACTER*1, INTENT(IN) :: mltype

REAL (wp), INTENT(OUT) :: vcov(nalpha,nalpha)

REAL (wp), INTENT(INOUT) :: log_lik

INTEGER, INTENT(INOUT) :: error_message

! Local variables

INTEGER :: i, ifail

REAL (wp) :: det, wk(nalpha+1,nalpha), wkspce(nalpha)

wk = 0.0_wp ! sum(t(x) *+ w *+ x)

DO i = 1, SIZE(haul)

wk(1:nalpha,:) = wk(1:nalpha,:) + MATMUL(TRANSPOSE(haul(i)%x), MATMUL(haul(i)%w, haul(i)%x))

END DO

ifail = 1 ! vcov = inverse(sum(t(x) *+ w *+ x))

CALL F01ABF(wk, nalpha+1, nalpha, vcov, nalpha, wkspce, ifail)

IF (ifail .ne. 0) THEN

CALL ERROR_CHECK(40+ifail, error_message)

RETURN

END IF

DO i = 1, nalpha ! fill in upper triangle

vcov(i,i:nalpha) = vcov(i:nalpha,i)

END DO

IF (mltype == 'r') THEN ! increment log-likelihood by log(det(sum(t(x) *+ w *+ x)))

CALL F03ABF(wk(1:nalpha,:), nalpha, nalpha, det, wkspce, ifail)

IF (ifail .ne. 0) THEN

CALL ERROR_CHECK(50+ifail, error_message)

RETURN

END IF

log_lik = log_lik + log(det)

END IF

END SUBROUTINE VCOV_CALC

! ****************************************************************************************************************************

SUBROUTINE B_CALC(haul, vcov, d, alpha, log_lik)

! calculate alpha, b, and increment the log-likelihood

IMPLICIT NONE

! Arguments

TYPE (HAUL_INFO), DIMENSION(:), INTENT(INOUT) :: haul

REAL (wp), INTENT(IN) :: vcov(nalpha,nalpha)

REAL (wp), INTENT(IN) :: d(nnu,nnu)

REAL (wp), INTENT(OUT) :: alpha(nalpha)

REAL (wp), INTENT(INOUT) :: log_lik

! Local variables

INTEGER :: i

REAL (wp) :: wk1(nnu), wk2(nnu)

alpha = 0.0_wp ! alpha = vcov *+ sum(t(x) *+ w *+ nu)

DO i = 1, SIZE(haul)

alpha = alpha + MATMUL(TRANSPOSE(haul(i)%x), MATMUL(haul(i)%w, haul(i)%nu))

END DO

alpha = MATMUL(vcov, alpha)

DO i = 1, SIZE(haul) ! b = d *+ w *+ (nu - x *+ alpha) and increment log-lik

wk1 = haul(i)%nu - MATMUL(haul(i)%x, alpha) ! by sum(t(nu - x *+ alpha) *+ w *+ (nu - x *+ alpha))

wk2 = MATMUL(haul(i)%w, wk1)

haul(i)%b = MATMUL(d, wk2)

log_lik = log_lik + SUM(wk1 * wk2)

END DO

END SUBROUTINE B_CALC

! ****************************************************************************************************************************

SUBROUTINE D_CALC(haul, vcov, mltype, d)

! update between-haul variance d

IMPLICIT NONE

! Arguments

TYPE (HAUL_INFO), DIMENSION(:), INTENT(IN) :: haul

REAL (wp), INTENT(IN) :: vcov(nalpha,nalpha)

CHARACTER*1, INTENT(IN) :: mltype

REAL (wp), INTENT(INOUT) :: d(nnu,nnu)

! Local variables

INTEGER :: nhaul, i

REAL (wp), DIMENSION(SIZE(haul), nnu) :: bwk

REAL (wp), DIMENSION(nnu,nnu) :: dnew, wk

! ML: calculate dnew = dold - sum(dold *+ w *+ d) / nhaul + sum(b *+ t(b)) / nhaul

nhaul = SIZE(haul)

wk = 0.0_wp

DO i = 1, nhaul

bwk(i,:) = haul(i)%b

wk = wk + haul(i)%w

END DO

wk = - MATMUL(d, MATMUL(wk, d)) + MATMUL(TRANSPOSE(bwk), bwk)

dnew = d + wk / nhaul

! REML: dnew = as above + sum(dold *+ w *+ x *+ vcov *+ t(x) *+ w *+ dold) / nhaul

IF (mltype .eq. 'r') THEN

wk = 0.0_wp

DO i = 1, nhaul

wk = wk + MATMUL(haul(i)%w, MATMUL(haul(i)%x, MATMUL(vcov, MATMUL(TRANSPOSE(haul(i)%x), haul(i)%w))))

END DO

dnew = dnew + MATMUL(d, MATMUL(wk, d)) / nhaul

END IF

d = dnew

RETURN

END SUBROUTINE D_CALC

! ****************************************************************************************************************************

SUBROUTINE MPRINT(log_lik, alpha, d)

! monitor convergence

IMPLICIT NONE

! Arguments

REAL(wp), INTENT(IN) :: log_lik

REAL(wp), INTENT(IN) :: alpha(nalpha)

REAL(wp), INTENT(IN) :: d(nnu,nnu)

! Local variables

INTEGER i, j

WRITE(*,'(//,A,//,A,e15.8)') '*** convergence monitoring ***', 'log-likelihood = ', log_lik

WRITE(*,'(/,A,/)') ' alpha'

DO i = 1, nalpha

WRITE(*,'(f10.4)') alpha(i)

END DO

WRITE(*,'(/,A,/)') ' d'

DO i = 1, nnu

WRITE(*,'(8f10.4)') (d(i,j), j = 1, nnu)

END DO

RETURN

END SUBROUTINE MPRINT

! ***************************************************************************************************************************

SUBROUTINE ERROR_CHECK(icheck, error_message)

IMPLICIT NONE

INTEGER, INTENT(IN) :: icheck

INTEGER, OPTIONAL, INTENT(OUT) :: error_message

IF (PRESENT(error_message)) THEN

error_message = icheck

RETURN

ELSE

SELECT CASE (icheck)

CASE (1)

WRITE(*,111) '*** ERROR ran out or iterations: error_message = ', icheck

CASE (2)

WRITE(*,111) '*** ERROR mltype is invalid: error_message = ', icheck

CASE (3)

WRITE(*,111) '*** ERROR d is not positive definite: error_message = ', icheck

CASE (4)

WRITE(*,111) '*** ERROR d is not specified: error_message = ', icheck

CASE (11:14)

WRITE(*,111) '*** ERROR calculating det(sum(t(x) *+ x)): NAG F03ABF ifail = ', icheck - 10

CASE (21:22)

WRITE(*,111) '*** ERROR calculating w matrix for some haul: NAG F01ABF ifail = ', icheck - 20

CASE (31:34)

WRITE(*,111) '*** ERROR calculating det(v) for some haul: NAG F03ABF ifail = ', icheck - 30

CASE (41:42)

WRITE(*,111) '*** ERROR calculating variance(alpha): NAG F01ABF ifail = ', icheck - 40

CASE (51:54)

WRITE(*,111) '*** ERROR calculating det(information matrix): NAG F03ABF ifail = ', icheck - 50

END SELECT

STOP

END IF

RETURN

111 FORMAT(//,A,i1,//)

END SUBROUTINE ERROR_CHECK

END SUBROUTINE FRSEL_FIT

END MODULE FRSEL