diff --git a/src/fstats_distributions.f90 b/src/fstats_distributions.f90 index bc30726..289cc46 100644 --- a/src/fstats_distributions.f90 +++ b/src/fstats_distributions.f90 @@ -31,6 +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. end type interface @@ -135,8 +137,40 @@ pure function distribution_property(this) result(rst) procedure, public :: variance => bd_variance end type -! ------------------------------------------------------------------------------ contains +! ------------------------------------------------------------------------------ +pure elemental function dist_area(this, x) result(rst) + !! Computes the area under the PDF curve up to the value of X specified. + class(distribution), intent(in) :: this + !! The distribution object. + real(real64), intent(in) :: x + !! The upper parameter limit. + real(real64) :: rst + !! The requested area. + + ! Local Variables + integer(int32), parameter :: maxiter = 100 + real(real64), parameter :: tol = 1.0d-6 + integer(int32) :: i + real(real64) :: f, df, h, twoh, dy + + ! Process + ! + ! We use a simplified Newton's method to solve for the independent variable + ! of the CDF function + h = 1.0d-6 + twoh = 2.0d0 * h + rst = 0.5d0 ! just an initial guess + do i = 1, maxiter + ! Compute the CDF and its derivative at y + f = this%cdf(rst) - x + df = (this%cdf(rst + h) - this%cdf(rst - h)) / twoh + dy = f / df + rst = rst - dy + if (abs(dy) < tol) exit + end do +end function + ! ****************************************************************************** ! STUDENT'S T-DISTRIBUTION ! ------------------------------------------------------------------------------ diff --git a/src/fstats_hypothesis.f90 b/src/fstats_hypothesis.f90 index cbfa664..9e3d86d 100644 --- a/src/fstats_hypothesis.f90 +++ b/src/fstats_hypothesis.f90 @@ -42,29 +42,11 @@ pure function confidence_interval_scalar(dist, alpha, s, n) result(rst) !! The result. ! Local Variables - integer(int32), parameter :: maxiter = 100 - real(real64), parameter :: tol = 1.0d-6 - integer(int32) :: i - real(real64) :: x, f, df, h, twoh, dy + real(real64) :: x ! Process - ! - ! We use a simplified Newton's method to solve for the independent variable - ! of the CDF function where it equals 1 - alpha / 2. - h = 1.0d-6 - twoh = 2.0d0 * h x = 1.0d0 - alpha / 2.0d0 - rst = 0.5d0 - do i = 1, maxiter - ! Compute the CDF and its derivative at y - f = dist%cdf(rst) - x - df = (dist%cdf(rst + h) - dist%cdf(rst - h)) / twoh - dy = f / df - rst = rst - dy - if (abs(dy) < tol) exit - end do - - ! Determine the actual interval + rst = dist%area(x) rst = rst * s / sqrt(real(n, real64)) end function @@ -468,5 +450,23 @@ subroutine levenes_test(x, stat, p, err) p = 1.0d0 - dist%cdf(stat) end subroutine +! ------------------------------------------------------------------------------ +pure function sample_size(dist, delta, pwr, 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. + class(distribution), intent(in) :: dist + !! The distribution to utilize as a measure. + 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 :: alpha + !! The desired significance level. The default for this value is 0.05. + integer(int32) :: rst + !! The minimum sample size requried to achieve the desired experimental + !! outcome. +end function + ! ------------------------------------------------------------------------------ end module \ No newline at end of file