s_lin_lda_2c Subroutine

public subroutine s_lin_lda_2c(x, nd, nv, nc, sa, g, score, mh)

2-class multivariate Linear Discriminant Analysis (LDA)

Performs classification and returns: - Standardised discriminant coefficients - Reclassification accuracy - Mahalanobis distance - Discriminant criterion

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: x(nd,nv,nc)

input data (nd samples × nv variables × nc classes)
integer(kind=i4), intent(in) :: nd

number of datapoints per class
integer(kind=i4), intent(in) :: nv

number of variables
integer(kind=i4), intent(in) :: nc

number of classes (must be 2)
real(kind=wp), intent(out) :: sa(nv)

standardised discriminant coefficients
real(kind=wp), intent(out) :: g

discriminant criterion
real(kind=wp), intent(out) :: score

classification score
real(kind=wp), intent(out), optional :: mh

Mahalanobis distance

Calls

proc~~s_lin_lda_2c~~CallsGraph proc~s_lin_lda_2c s_lin_lda_2c eigh eigh proc~s_lin_lda_2c->eigh proc~f_sts_cov_core f_sts_cov_core proc~s_lin_lda_2c->proc~f_sts_cov_core proc~f_sts_mean_core f_sts_mean_core proc~s_lin_lda_2c->proc~f_sts_mean_core proc~f_sts_var_core f_sts_var_core proc~s_lin_lda_2c->proc~f_sts_var_core proc~s_err_print s_err_print proc~s_lin_lda_2c->proc~s_err_print proc~f_sts_cov_core->proc~f_sts_mean_core proc~f_sts_var_core->proc~f_sts_mean_core proc~f_utl_r2c f_utl_r2c proc~s_err_print->proc~f_utl_r2c

Called by

proc~~s_lin_lda_2c~~CalledByGraph proc~s_lin_lda_2c s_lin_lda_2c interface~fsml_lda_2class fsml_lda_2class interface~fsml_lda_2class->proc~s_lin_lda_2c

Source Code

subroutine s_lin_lda_2c(x, nd, nv, nc, sa, g, score, mh)

! ==== Description
!! 2-class multivariate Linear Discriminant Analysis (LDA)
!!
!! Performs classification and returns:
!! - Standardised discriminant coefficients
!! - Reclassification accuracy
!! - Mahalanobis distance
!! - Discriminant criterion

! ==== Declarations
  integer(i4), intent(in)            :: nc              !! number of classes (must be 2)
  integer(i4), intent(in)            :: nv              !! number of variables
  integer(i4), intent(in)            :: nd              !! number of datapoints per class
  real(wp)   , intent(in)            :: x(nd,nv,nc)     !! input data (nd samples × nv variables × nc classes)
  real(wp)   , intent(out)           :: score           !! classification score
  real(wp)   , intent(out)           :: sa(nv)          !! standardised discriminant coefficients
  real(wp)   , intent(out)           :: g               !! discriminant criterion
  real(wp)   , intent(out), optional :: mh              !! Mahalanobis distance
  real(wp)                           :: xmv(nv,nc)      !! group mean vectors
  real(wp)                           :: s_g(nv,nv,nc)   !! group covariance matrices
  real(wp)                           :: s_pool(nv,nv)   !! pooled covariance matrix
  real(wp)                           :: s_pool_i(nv,nv) !! inverse of pooled covariance matrix
  real(wp)                           :: a(nv)           !! discriminant vector
  real(wp)                           :: d_pool(nc*nd)   !! pooled data for std calc
  real(wp)                           :: tmpv(nd)        !! temporary vector
  real(wp)                           :: tmp             !! temporary scalars
  integer(i4)                        :: i, j, k         !! loop counters
  real(wp)                           :: ew(nv)          !! eigenvalues of pooled covariance
  real(wp)                           :: ev(nv,nv)       !! eigenvectors of pooled covariance
  real(wp)                           :: ew_diag(nv,nv)  !! diagonal matrix of inverted eigenvalues

! ==== Instructions

  ! ---- validate inputs
  if (nc .ne. 2) then
     call s_err_print("[fsml error] 2-class LDA: Number of classes must be 2.")
     error stop
  endif
  if (nv .lt. 2) then
     call s_err_print("[fsml error] 2-class LDA: 2+ variables required.")
     error stop
  endif

  ! ---- compute group means and covariance matrices
  xmv = 0.0_wp
  s_g = 0.0_wp
  do i = 1, nc
     do j = 1, nv
        tmpv(:) = x(:,j,i)
        xmv(j,i) = f_sts_mean_core(tmpv)
     enddo
     do j = 1, nv
        do k = 1, nv
           ! get sample covariance (ddf set to 1)
           s_g(k,j,i) = f_sts_cov_core(x(:,j,i), x(:,k,i), ddf=1.0_wp)
        enddo
     enddo
  enddo

  ! ---- invert pooled covariance matrix using eigen-decomposition

  ! compute pooled covariance matrix (equal nd → unweighted average)
  s_pool = 0.5_wp * (s_g(:,:,1) + s_g(:,:,2))

  ! eigendecomposition of s_pool
  call eigh(s_pool, ew, vectors=ev)

  ! set returns to sentinel
  do i = 1, nv
     if (ew(i) .le. 0.0_wp) then
        if (present(mh)) mh = c_sentinel_r
        g     = c_sentinel_r
        sa(:) = c_sentinel_r
        score = c_sentinel_r
        return
     endif
  enddo

  ! construct diagonal matrix from inverted eigenvalues
  ew_diag = 0.0_wp
  do i = 1, nv
     ew_diag(i,i) = 1.0_wp / ew(i)
  enddo

  ! compute inverse of pooled covariance matrix: V * D⁻¹ * Vᵗ
  s_pool_i = matmul(ev, matmul(ew_diag, transpose(ev)))

  ! ---- discriminant coefficients

  ! compute discriminant vector a = S_pool⁻¹ * (μ1 - μ2)
  a = matmul(s_pool_i, xmv(:,1) - xmv(:,2))

  ! standardise coefficients
  do i = 1, nv
     d_pool(1:nd)       = x(1:nd,i,1)
     d_pool(nd+1:nc*nd) = x(1:nd,i,2)
     ! get sample variance (ddf set to 1)
     tmp   = f_sts_var_core(d_pool, ddf=1.0_wp)
     sa(i) = a(i) * sqrt(tmp)
  enddo

  ! ---- compute Mahalanobis distance
  if (present(mh)) then
     mh = sqrt( dot_product( xmv(:,1) - xmv(:,2), &
          & matmul( s_pool_i, xmv(:,1) - xmv(:,2) ) ) )
  endif

  ! ---- compute discriminant criterion g
  g = 0.5_wp * (dot_product(a, xmv(:,1)) + dot_product(a, xmv(:,2)))

  ! ---- re-classification and scoring
  score = 0.0_wp
  do i = 1, nd
     do j = 1, nc
        tmp = dot_product(a, x(i,:,j))
        if ((tmp .ge. g .and. j .eq. 1) .or. &
         & (tmp .lt. g .and. j .eq. 2)) then
           score = score + 1.0_wp
        endif
     enddo
  enddo
  score = score / real(nc * nd, kind=wp)

end subroutine s_lin_lda_2c