Skip to content

Commit

Permalink
The commit with sha cbe0ca on 5 April 2023 transitioned to the use of…
Browse files Browse the repository at this point in the history 
… distribution_param_types in

the probit_transform_mod.f90 to communicate parameters of arbitrary distributions to the individual
distribution modules. This was required for and tested in the bnrh_distribution_mod.f90. However,
it also had a footprint in the beta_distribution_mod.f90 and the gamma_distribution_mod.f90 and
these changes were not tested. In addition, this commit provided partial support for a generalized
beta distribution in which the bounds are not 0 and 1. It also provided partial support for a
generalized gamma distribution where the lower bound did not necessarily have to be 0 and where the
gamma could be bounded above rather than below. Neither of these generalizations was ever completed.
The use of the standard test subroutine in the beta_distribution_mod would have revealed that there
was no longer a guarantee that the appropriate parameters were set for the beta module and one of
the tests would have failed; that test still fails in the commit prior to this one. The gamma
distribution had a similar problem, all parameters were not set, but the result does not lead to a
failure of the gamma test subroutine.

This commit fixes the errors in the beta_distribution_mod.f90 and gamma_distribution_mod.f90. It
also upgrades the test subroutines in these modules and the normal_distribution_mod.f90 to comply
with DART testing standards. These subroutines now explicitly print ‘PASS’ for tests that pass and
‘FAIL’ for tests that fail allowing more automated testing scripts. It also creates directories in
the developer_tests directory with programs that run the test subroutine for each of the
distributions. Details of differences in each module follow.

_________________________________________________________________
Beta_distribution_mod.f90: The key change is to remove partial support for generalized beta
distributions. A comment now makes it clear that this module only supports the standard beta
distribution that is bounded by 0 and 1. Subroutines beta_cdf, inverse_beta_cdf,
set_beta_params_from_ens, and inv_beta_first_guess no longer have the lower_bound and upper_bound
as arguments (they are only allowed to be 0 and 1). In addition, inv_beta_first_guess no longer has
bounded_below and bounded_above as arguments as they must both be true. inv_beta_cdf now sets the
bounded_below, bounded_above, lower_bound and upper_bound parameters to the correct values for the
beta. The partially implemented distribution scaling was removed from inv_beta_cdf.

Functions inv_beta_cdf_params and beta_cdf_params now include an error check to make sure that the
lower and upper bounds in the distribution_params_type have been set to 0 and 1 as other values are
not supported.

Subroutine set_beta_params_from_ens has changed the distribution_params_type to an intent out
argument and defines all six parameters correctly. It also set the parameter type to
BETA_DISTRIBUTION.

Subroutine test_beta has been modified to print the mandated ‘PASS’ or ‘FAIL’ as appropriate.

_________________________________________________________________
gamma_distribution_mod.f90: The key change is to remove partial support for generalized gamma
distributions. A comment now makes it clear that this module only supports the standard gamma
distribution that is bounded below by 0 and unbounded above. Subroutines gamma_cdf, inverse_gamma_cdf,
set_gamma_params_from_ens, and inv_gamma_first_guess no longer have bounded_below, bounded_above,
lower_bound, and upper_bound as arguments. inv_gamma_cdf now sets the bounded_below, bounded_above,
lower_bound and upper_bound parameters to the correct values for the gamma.

Functions inv_gamma_cdf_params and gamma_cdf_params now include an error check to make sure that the
lower and upper bounds in the distribution_params_type have been set to 0 and missing_r8 as other
values are not supported.

Subroutine set_gamma_params_from_ens has changed the distribution_params_type to an intent out
argument and defines all six parameters correctly. It also sets the parameter type to
GAMMA_DISTRIBUTION.

Subroutine test_gamma has been modified to print the mandated ‘PASS’ or ‘FAIL’ as appropriate.

_________________________________________________________________
normal_distribution_mod.f90: There were no known errors in this module. Several changes were made
to be consistent with the appropriate use of the distribution_params_type.

Functions inv_std_normal_cdf and set_normal_params_from_ens now set the appropriate values for the
bounded_below, bounded_above, lower_bound and upper_bound components of the distribution_params_type.
The distribution_params_type was changed to intent out in set_normal_params_from_ens.

The magic number definition of the maximum delta in the inv_cdf root searching routine was changed
to be a parameter but the value of 1e-8 was unchanged. A comment notes that changing this parameter to
1e-9 allows all of Ian Groom’s KDE distribution tests to pass.

Subroutine test_normal was modified to print the mandated ‘PASS’ or ‘FAIL’ as appropriate.
_________________________________________________________________
assim_tools_mod.f90: Deprecated arguments in calls to gamma_cdf and inv_gamma_cdf were removed. Note
that the beta distribution has not been implemented as a prior for observation space in
assim_tools_mod.f90 so there is no corresponding change.

_________________________________________________________________
probit_transform_mod.f90: Deprecated arguments were removed from a call to set_gamma_params_from_ens
and set_beta_params_from_ens. Arguments lower_bound and upper_bound were removed from subroutine
to_probit_beta.

_________________________________________________________________
Directories normal_distribution, gamma_distribution, and beta_distribution were created in
developer_tests. They contain programs test_normal_dist.f90, test_gamma_dist.f90, and
test_beta_dist.f90 and a work subdirectory. The work subdirectories contain a quickbuild.sh and an
input.nml. Each program runs the corresponding test_ subroutine for its distribution. All tests are
PASS with this commit.
  • Loading branch information
@jlaucar
jlaucar committed Aug 20, 2024
1 parent f193ab4 commit 5b3850f
Show file tree
Hide file tree
Showing 14 changed files with 339 additions and 81 deletions.
4 changes: 2 additions & 2 deletions assimilation_code/modules/assimilation/assim_tools_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1057,7 +1057,7 @@ subroutine obs_increment_gamma(ens, ens_size, prior_mean, prior_var, obs, obs_va
! Compute the prior quantiles of each ensemble member in the prior gamma distribution
call gamma_mn_var_to_shape_scale(prior_mean, prior_var, prior_shape, prior_scale)
do i = 1, ens_size
q(i) = gamma_cdf(ens(i), prior_shape, prior_scale, .true., .false., 0.0_r8, missing_r8)
q(i) = gamma_cdf(ens(i), prior_shape, prior_scale)
end do

! Compute the statistics of the continous posterior distribution
Expand All @@ -1074,7 +1074,7 @@ subroutine obs_increment_gamma(ens, ens_size, prior_mean, prior_var, obs, obs_va

! Now invert the quantiles with the posterior distribution
do i = 1, ens_size
post(i) = inv_gamma_cdf(q(i), post_shape, post_scale, .true., .false., 0.0_r8, missing_r8)
post(i) = inv_gamma_cdf(q(i), post_shape, post_scale)
end do

obs_inc = post - ens
Expand Down
94 changes: 61 additions & 33 deletions assimilation_code/modules/assimilation/beta_distribution_mod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,9 @@

! Thanks to Chris Riedel who developed the methods in this module.

! This module only supports standard beta distributions for which the
! lower bound is 0 and the upper bound is 1.

module beta_distribution_mod

use types_mod, only : r8, PI, missing_r8
Expand All @@ -12,7 +15,7 @@ module beta_distribution_mod

use random_seq_mod, only : random_seq_type, random_uniform

use distribution_params_mod, only : distribution_params_type
use distribution_params_mod, only : distribution_params_type, BETA_DISTRIBUTION

use normal_distribution_mod, only : inv_cdf

Expand All @@ -36,11 +39,11 @@ subroutine test_beta

! This routine provides limited tests of the numerics in this module. It begins
! by comparing a handful of cases of the pdf and cdf to results from Matlab. It
! then tests the quality of the inverse cdf for a single shape/scale pair. Failing
! then tests the quality of the inverse cdf for a single alpha/beta pair. Failing
! these tests suggests a serious problem. Passing them does not indicate that
! there are acceptable results for all possible inputs.

real(r8) :: x, y, p, inv
real(r8) :: x, y, inv
real(r8) :: alpha, beta, max_diff
integer :: i

Expand All @@ -62,9 +65,14 @@ subroutine test_beta
write(*, *) 'Absolute value of differences should be less than 1e-15'
do i = 1, 7
pdf_diff(i) = beta_pdf(mx(i), malpha(i), mbeta(i)) - mpdf(i)
cdf_diff(i) = beta_cdf(mx(i), malpha(i), mbeta(i), 0.0_r8, 1.0_r8) - mcdf(i)
cdf_diff(i) = beta_cdf(mx(i), malpha(i), mbeta(i)) - mcdf(i)
write(*, *) i, pdf_diff(i), cdf_diff(i)
end do
if(maxval(pdf_diff) < 1e-15_r8 .and. maxval(cdf_diff) < 1e-15_r8) then
write(*, *) 'Matlab Comparison Tests: PASS'
else
write(*, *) 'Matlab Comparison Tests: FAIL'
endif

! Test many x values for cdf and inverse cdf for a single set of alpha and beta
alpha = 5.0_r8
Expand All @@ -73,15 +81,19 @@ subroutine test_beta
max_diff = -1.0_r8
do i = 0, 1000
x = i / 1000.0_r8
p = beta_pdf(x, alpha, beta)
y = beta_cdf(x, alpha, beta, 0.0_r8, 1.0_r8)
inv = inv_beta_cdf(y, alpha, beta, 0.0_r8, 1.0_r8)
y = beta_cdf(x, alpha, beta)
inv = inv_beta_cdf(y, alpha, beta)
max_diff = max(abs(x - inv), max_diff)
end do

write(*, *) '----------------------------'
write(*, *) 'max difference in inversion is ', max_diff
write(*, *) 'max difference should be less than 1e-14'
if(max_diff < 1e-14_r8) then
write(*, *) 'Inversion Tests: PASS'
else
write(*, *) 'Inversion Tests: FAIL'
endif

end subroutine test_beta

Expand All @@ -93,20 +105,25 @@ function inv_beta_cdf_params(quantile, p) result(x)
real(r8), intent(in) :: quantile
type(distribution_params_type), intent(in) :: p

! Only standard beta is currently supported. Fail if bounds are not 0 and 1
if(p%lower_bound /= 0.0_r8 .or. p%upper_bound /= 1.0_r8) then
errstring = 'Only standard beta distribution with bounds of 0 and 1 is supported'
call error_handler(E_ERR, 'inv_beta_cdf_params', errstring, source)
endif

x = inv_cdf(quantile, beta_cdf_params, inv_beta_first_guess_params, p)

end function inv_beta_cdf_params

!-----------------------------------------------------------------------

function inv_beta_cdf(quantile, alpha, beta, lower_bound, upper_bound) result(x)
function inv_beta_cdf(quantile, alpha, beta) result(x)

real(r8) :: x
real(r8), intent(in) :: quantile
real(r8), intent(in) :: alpha, beta
real(r8), intent(in) :: lower_bound, upper_bound

! Given a quantile, finds the value of x for which the scaled beta cdf
! Given a quantile, finds the value of x for which the beta cdf
! with alpha and beta has approximately this quantile

type(distribution_params_type) :: p
Expand All @@ -116,15 +133,13 @@ function inv_beta_cdf(quantile, alpha, beta, lower_bound, upper_bound) result(x)
call error_handler(E_ERR, 'inv_beta_cdf', errstring, source)
endif

! Load the p type for the generic cdf calls
p%params(1) = alpha; p%params(2) = beta
! Beta must be bounded on both sides
p%lower_bound = lower_bound; p%upper_bound = upper_bound
p%bounded_below = .true.; p%bounded_above = .true.
p%lower_bound = 0.0_r8; p%upper_bound = 1.0_r8

x = inv_beta_cdf_params(quantile, p)

! Undo the scaling
x = x * (upper_bound - lower_bound) + lower_bound

end function inv_beta_cdf

!---------------------------------------------------------------------------
Expand Down Expand Up @@ -171,25 +186,34 @@ function beta_cdf_params(x, p)
real(r8), intent(in) :: x
type(distribution_params_type), intent(in) :: p

! A translation routine that is required to use the generic cdf optimization routine
! Extracts the appropriate information from the distribution_params_type that is needed
! for a call to the function beta_cdf below.

real(r8) :: alpha, beta

! Only standard beta is currently supported. Fail if bounds are not 0 and 1
if(p%lower_bound /= 0.0_r8 .or. p%upper_bound /= 1.0_r8) then
errstring = 'Only standard beta distribution with bounds of 0 and 1 is supported'
call error_handler(E_ERR, 'beta_cdf_params', errstring, source)
endif

alpha = p%params(1); beta = p%params(2)
beta_cdf_params = beta_cdf(x, alpha, beta, p%lower_bound, p%upper_bound)
beta_cdf_params = beta_cdf(x, alpha, beta)

end function beta_cdf_params

!---------------------------------------------------------------------------

function beta_cdf(x, alpha, beta, lower_bound, upper_bound)
function beta_cdf(x, alpha, beta)

! Returns the cumulative distribution of a beta function with alpha and beta
! at the value x

! Returns a large negative value if called with illegal values
! Returns a failed_value if called with illegal values

real(r8) :: beta_cdf
real(r8), intent(in) :: x, alpha, beta
real(r8), intent(in) :: lower_bound, upper_bound

! Parameters must be positive
if(alpha <= 0.0_r8 .or. beta <= 0.0_r8) then
Expand Down Expand Up @@ -238,7 +262,7 @@ function random_beta(r, alpha, beta)
! Draw from U(0, 1) to get a quantile
quantile = random_uniform(r)
! Invert cdf to get a draw from beta
random_beta = inv_beta_cdf(quantile, alpha, beta, 0.0_r8, 1.0_r8)
random_beta = inv_beta_cdf(quantile, alpha, beta)

end function random_beta

Expand Down Expand Up @@ -326,24 +350,25 @@ function inv_beta_first_guess_params(quantile, p)
real(r8), intent(in) :: quantile
type(distribution_params_type), intent(in) :: p

! A translation routine that is required to use the generic first_guess for
! the cdf optimization routine.
! Extracts the appropriate information from the distribution_params_type that is needed
! for a call to the function inv_beta_first_guess below.

real(r8) :: alpha, beta

alpha = p%params(1); beta = p%params(2)
inv_beta_first_guess_params = inv_beta_first_guess(quantile, alpha, beta, &
p%bounded_below, p%bounded_above, p%lower_bound, p%upper_bound)
inv_beta_first_guess_params = inv_beta_first_guess(quantile, alpha, beta)

end function inv_beta_first_guess_params

!---------------------------------------------------------------------------

function inv_beta_first_guess(x, alpha, beta, &
bounded_below, bounded_above, lower_bound, upper_bound)
function inv_beta_first_guess(x, alpha, beta)

real(r8) :: inv_beta_first_guess
real(r8), intent(in) :: x
real(r8), intent(in) :: alpha, beta
logical, intent(in) :: bounded_below, bounded_above
real(r8), intent(in) :: lower_bound, upper_bound

! Need some sort of first guess, should be smarter here
! For starters, take the mean for this alpha and beta
Expand Down Expand Up @@ -376,19 +401,22 @@ end subroutine beta_alpha_beta

!---------------------------------------------------------------------------

subroutine set_beta_params_from_ens(ens, num, lower_bound, upper_bound, p)
subroutine set_beta_params_from_ens(ens, num, p)

integer, intent(in) :: num
real(r8), intent(in) :: ens(num)
real(r8), intent(in) :: lower_bound, upper_bound
type(distribution_params_type), intent(inout) :: p
integer, intent(in) :: num
real(r8), intent(in) :: ens(num)
type(distribution_params_type), intent(out) :: p

real(r8) :: alpha, beta

! Set up the description of the beta distribution defined by the ensemble
p%distribution_type = BETA_DISTRIBUTION

! Set the bounds info
p%lower_bound = lower_bound; p%upper_bound = upper_bound
p%bounded_below = .true.; p%bounded_above = .true.
p%lower_bound = 0.0_r8; p%upper_bound = 1.0_r8

! Get alpha and beta for the scaled ensemble
! Get alpha and beta for the ensemble
call beta_alpha_beta(ens, num, alpha, beta)
p%params(1) = alpha
p%params(2) = beta
Expand Down
Loading

0 comments on commit 5b3850f

Please sign in to comment.