Add trimmed mean calculation

src/fstats.f90

@@ -49,6 +49,7 @@ module fstats
     public :: jacobian
     public :: nonlinear_least_squares
     public :: allan_variance
+    public :: trimmed_mean
     public :: FS_LEVENBERG_MARQUARDT_UPDATE
     public :: FS_QUADRATIC_UPDATE
     public :: FS_NIELSEN_UPDATE
@@ -691,6 +692,12 @@ pure module function cs_variance(this) result(rst)
         module procedure :: confidence_interval_real32_array
     end interface
 
+    interface trimmed_mean
+        !! Computes the trimmed mean of a data set.
+        module procedure :: trimmed_mean_real64
+        module procedure :: trimmed_mean_real32
+    end interface
+
     ! statistics_implementation.f90
     interface
         pure module function mean_real64(x) result(rst)
@@ -1089,6 +1096,32 @@ module function anova_model_fit(nmodelparams, ymeas, ymod, err) &
                 !! - FS_MEMORY_ERROR: Occurs if a memory error is encountered.
             type(single_factor_anova_table) :: rst
         end function
+
+        module function trimmed_mean_real64(x, p) result(rst)
+            !! Computes the trimmed mean of a data set.
+            real(real64), intent(inout), dimension(:) :: x
+                !! An N-element array containing the data.  On output, the
+                !! array is sorted into ascending order.
+            real(real64), intent(in), optional :: p
+                !! An optional parameter specifying the percentage of values
+                !! from either end of the distribution to remove.  The default
+                !! is 0.05 such that the bottom 5% and top 5% are removed.
+            real(real64) :: rst
+                !! The trimmed mean.
+        end function
+
+        module function trimmed_mean_real32(x, p) result(rst)
+            !! Computes the trimmed mean of a data set.
+            real(real32), intent(inout), dimension(:) :: x
+                !! An N-element array containing the data.  On output, the
+                !! array is sorted into ascending order.
+            real(real32), intent(in), optional :: p
+                !! An optional parameter specifying the percentage of values
+                !! from either end of the distribution to remove.  The default
+                !! is 0.05 such that the bottom 5% and top 5% are removed.
+            real(real32) :: rst
+                !! The trimmed mean.
+        end function
     end interface
 
 ! ******************************************************************************

src/statistics_implementation.f90

@@ -1057,10 +1057,66 @@ module function anova_model_fit(nmodelparams, ymeas, ymod, err) result(rst)
 101 format(A, I0, A)
 end function
 
+! ------------------------------------------------------------------------------
+module function trimmed_mean_real64(x, p) result(rst)
+    ! Arguments
+    real(real64), intent(inout), dimension(:) :: x
+    real(real64), intent(in), optional :: p
+    real(real64) :: rst
+
+    ! Local Variables
+    integer(int32) :: i1, i2, n
+    real(real64) :: pv
+
+    ! Initialization
+    if (present(p)) then
+        pv = abs(p)
+    else
+        pv = 0.05d0
+    end if
+
+    ! Sort the array into ascending order
+    call sort(x, .true.)
+
+    ! Find the limiting indices
+    n = size(x)
+    i1 = max(floor(n * pv, int32), 1)
+    i2 = min(n, n - i1 + 1)
+    rst = mean(x(i1:i2))
+end function
+
+! --------------------
+module function trimmed_mean_real32(x, p) result(rst)
+    ! Arguments
+    real(real32), intent(inout), dimension(:) :: x
+    real(real32), intent(in), optional :: p
+    real(real32) :: rst
+
+    ! Local Variables
+    integer(int32) :: i1, i2, n
+    real(real32) :: pv
+
+    ! Initialization
+    if (present(p)) then
+        pv = abs(p)
+    else
+        pv = 0.05
+    end if
+
+    ! Sort the array into ascending order
+    call r32_sort(x)
+
+    ! Find the limiting indices
+    n = size(x)
+    i1 = max(floor(n * pv, int32), 1)
+    i2 = min(n, n - i1 + 1)
+    rst = mean(x(i1:i2))
+end function
+
 ! ******************************************************************************
 ! PRIVATE ROUTINES
 ! ------------------------------------------------------------------------------
-! This is a hack as the linalg library doesn't support 32-bit real32ing point
+! This is a hack as the linalg library doesn't support 32-bit floating point
 ! routines
 subroutine r32_sort(x)
     ! Arguments

tests/fstats_statistics_tests.f90

@@ -1014,5 +1014,45 @@ function incomplete_gamma_test_1() result(rst)
         end if
     end function
 
+! ------------------------------------------------------------------------------
+    function trimmed_mean_test_1() result(rst)
+        use linalg, only : sort
+
+        ! Arguments
+        logical :: rst
+
+        ! Parameters
+        integer(int32), parameter :: n = 100
+        integer(int32), parameter :: i1 = 5
+        integer(int32), parameter :: i2 = n - i1 + 1
+        real(real64), parameter :: p = 0.05d0
+
+        ! Local Variables
+        real(real64) :: x(n), tm, ans
+        real(real32) :: x32(n), tm32, ans32
+
+        ! Initialization
+        rst = .true.
+        call random_number(x)
+        call sort(x, .true.) ! sort into ascending order
+        x32 = real(x, real32)
+        ans = mean(x(i1:i2))
+        ans32 = mean(x32(i1:i2))
+
+        ! Test 1
+        tm = trimmed_mean(x, p = p)
+        if (.not.is_equal(tm, ans)) then
+            rst = .false.
+            print '(A)', "TEST FAILED: Trimmed Mean Test 1"
+        end if
+
+        ! Test 2
+        tm32 = trimmed_mean(x32, p = real(p, real32))
+        if (.not.is_equal(tm32, ans32)) then
+            rst = .false.
+            print '(A)', "TEST FAILED: Trimmed Mean Test 2"
+        end if
+    end function
+
 ! ------------------------------------------------------------------------------
 end module

tests/fstats_tests.f90

@@ -97,6 +97,10 @@ program tests
     local = test_allan_variance()
     if (.not.local) overall = .false.
 
+    ! Trimmed Mean Tests
+    local = trimmed_mean_test_1()
+    if (.not.local) overall = .false.
+
     ! End
     if (.not.overall) then
         stop 1