Update power calculations

src/fstats.f90

@@ -81,6 +81,7 @@ module fstats
     public :: pooled_variance
     public :: bartletts_test
     public :: levenes_test
+    public :: sample_size
     public :: FS_LEVENBERG_MARQUARDT_UPDATE
     public :: FS_QUADRATIC_UPDATE
     public :: FS_NIELSEN_UPDATE

src/fstats_distributions.f90

@@ -31,8 +31,8 @@ module fstats_distributions
             !! Computes the mode of the distribution.
         procedure(distribution_property), deferred, pass :: variance
             !! Computes the variance of the distribution.
-        procedure, public :: area => dist_area
-            !! Computes the area under the PDF curve up to the value specified.
+        procedure, public :: standardized_variable => dist_std_var
+            !! Computes the standardized variable for the distribution.
     end type
 
     interface
@@ -139,14 +139,14 @@ pure function distribution_property(this) result(rst)
 
 contains
 ! ------------------------------------------------------------------------------
-pure elemental function dist_area(this, x) result(rst)
-    !! Computes the area under the PDF curve up to the value of X specified.
+pure elemental function dist_std_var(this, x) result(rst)
+    !! Computes the standardized variable for the distribution.
     class(distribution), intent(in) :: this
         !! The distribution object.
     real(real64), intent(in) :: x
-        !! The upper parameter limit.
+        !! The value of interest.
     real(real64) :: rst
-        !! The requested area.
+        !! The result.
 
     ! Local Variables
     integer(int32), parameter :: maxiter = 100

src/fstats_hypothesis.f90

@@ -14,6 +14,7 @@ module fstats_hypothesis
     public :: f_test
     public :: bartletts_test
     public :: levenes_test
+    public :: sample_size
 
     interface confidence_interval
         !! Computes the confidence interval for the specified distribution.
@@ -46,7 +47,7 @@ pure function confidence_interval_scalar(dist, alpha, s, n) result(rst)
 
     ! Process
     x = 1.0d0 - alpha / 2.0d0
-    rst = dist%area(x)
+    rst = dist%standardized_variable(x)
     rst = rst * s / sqrt(real(n, real64))
 end function
 
@@ -451,21 +452,48 @@ subroutine levenes_test(x, stat, p, err)
 end subroutine
 
 ! ------------------------------------------------------------------------------
-pure function sample_size(dist, delta, pwr, alpha) result(rst)
+pure function sample_size(dist, var, delta, bet, alpha) result(rst)
     !! Estimates the sample size required to achieve an experiment with the
     !! desired power and significance levels to ascertain the desired 
     !! difference in parameter.
+    !!
+    !! See Also
+    !!
+    !! - [Wikipedia](https://en.wikipedia.org/wiki/Power_of_a_test)
     class(distribution), intent(in) :: dist
         !! The distribution to utilize as a measure.
+    real(real64), intent(in) :: var
+        !! An estimate of the population variance.
     real(real64), intent(in) :: delta
         !! The parameter difference that is desired.
-    real(real64), intent(in), optional :: pwr
-        !! The desired power level.  The default for this value is 0.8.
+    real(real64), intent(in), optional :: bet
+        !! The desired power level.  The default for this value is 0.2, for a 
+        !! power of 80%.
     real(real64), intent(in), optional :: alpha
-        !! The desired significance level.  The default for this value is 0.05.
-    integer(int32) :: rst
+        !! The desired significance level.  The default for this value is 0.05
+        !! for a confidence level of 95%.
+    real(real64) :: rst
         !! The minimum sample size requried to achieve the desired experimental
         !! outcome.
+
+    ! Local Variables
+    real(real64) :: a, b, za, zb
+
+    ! Initialization
+    if (present(bet)) then
+        b = bet
+    else
+        b = 0.8d0
+    end if
+    if (present(alpha)) then
+        a = alpha
+    else
+        a = 0.05d0
+    end if
+
+    za = dist%standardized_variable(1.0d0 - a / 2.0d0)
+    zb = dist%standardized_variable(b)
+    rst = 2.0d0 * (za + zb)**2 * var / (delta**2)
 end function
 
 ! ------------------------------------------------------------------------------

tests/fstats_distribution_tests.f90

@@ -186,5 +186,38 @@ function binomial_distribution_test_1() result(rst)
         end if
     end function
 
+! ------------------------------------------------------------------------------
+    function test_standardized_variable() result(rst)
+        ! Arguments
+        logical :: rst
+
+        ! Variables
+        real(real64), parameter :: tol = 1.0d-2
+        real(real64), parameter :: alpha1 = 0.975d0
+        real(real64), parameter :: alpha2 = 0.2d0
+        real(real64), parameter :: ans1 = 1.96d0
+        real(real64), parameter :: ans2 = 0.84d0
+        real(real64) :: z1, z2
+        type(normal_distribution) :: dist
+
+        ! Initialization
+        rst = .true.
+        call dist%standardize()
+
+        ! Test 1
+        z1 = dist%standardized_variable(alpha1)
+        if (.not.is_equal(z1, ans1, tol)) then
+            rst = .false.
+            print '(A)', "TEST FAILED: Standardized variable test -1."
+        end if
+
+        ! Test 2
+        z2 = dist%standardized_variable(0.8d0)
+        if (.not.is_equal(z2, ans2, tol)) then
+            rst = .false.
+            print '(A)', "TEST FAILED: Standardized variable test -2."
+        end if
+    end function
+
 ! ------------------------------------------------------------------------------
 end module

tests/fstats_statistics_tests.f90

@@ -1040,5 +1040,32 @@ function test_levene_1() result(rst)
         end if
     end function
 
+! ------------------------------------------------------------------------------
+    function test_sample_size() result(rst)
+        ! Arguments
+        logical :: rst
+
+        ! Variables
+        real(real64), parameter :: tol = 1.0d-3
+        real(real64), parameter :: ans = 15.698d0
+        real(real64), parameter :: var = 1.0d0
+        real(real64), parameter :: delta = 1.0d0
+        real(real64), parameter :: alpha = 0.05d0
+        real(real64), parameter :: bet = 0.8d0
+        real(real64) :: dn
+        type(normal_distribution) :: dist
+
+        ! Initialization
+        rst = .true.
+        call dist%standardize()
+
+        ! Test 1
+        dn = sample_size(dist, var, delta, bet, alpha)
+        if (.not.is_equal(ans, dn, tol)) then
+            rst = .false.
+            print '(A)', "TEST FAILED: Sample size test 1."
+        end if
+    end function
+
 ! ------------------------------------------------------------------------------
 end module

tests/fstats_tests.f90

@@ -128,6 +128,12 @@ program tests
     local = test_levene_1()
     if (.not.local) overall = .false.
 
+    local = test_standardized_variable()
+    if (.not.local) overall = .false.
+
+    local = test_sample_size()
+    if (.not.local) overall = .false.
+
     ! End
     if (.not.overall) then
         stop 1