Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Distribution fixes #788

Open
wants to merge 10 commits into
base: main
Choose a base branch
from
Open

Distribution fixes #788

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
9992f22
fix: beta_distrtibution_mod only supports standard beat distribution 0,1
hkershaw-brown Dec 26, 2024
be7f568
gamma_distribution_mod.f90: remove partial support for generalized ga…
hkershaw-brown Dec 26, 2024
b9e67b4
normal_distribution_mod: set distribution_params_types
hkershaw-brown Dec 26, 2024
86c7266
fix: any task with a failure to converge prints a message
hkershaw-brown Dec 26, 2024
55b7042
have UNSET as the default initalization for distribution_params_type
hkershaw-brown Jan 6, 2025
538d238
doc: note on standard bounds only for GAMMA and BETA
hkershaw-brown Jan 7, 2025
e0fb360
feat: force standard gamma and beta in distribution in qceff table read
hkershaw-brown Jan 7, 2025
79545d6
Merge branch 'main' into distribution-fixes
hkershaw-brown Jan 9, 2025
7930bff
set the distribution type for normal distribution
hkershaw-brown Jan 9, 2025
bcbeefe
remove unused bound arguments for to_probit_gamma
hkershaw-brown Jan 9, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions assimilation_code/modules/assimilation/assim_tools_mod.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1065,7 +1065,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 @@ -1082,7 +1082,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
91 changes: 53 additions & 38 deletions assimilation_code/modules/assimilation/beta_distribution_mod.f90
View file
Open in desktop
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 @@ -61,41 +64,37 @@ subroutine test_beta
! Compare to matlab
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)
pdf_diff(i) = abs(beta_pdf(mx(i), malpha(i), mbeta(i)) - mpdf(i))
cdf_diff(i) = abs(beta_cdf(mx(i), malpha(i), mbeta(i)) - mcdf(i))
write(*, *) i, pdf_diff(i), cdf_diff(i)
end do
if(abs(maxval(pdf_diff)) < 1e-15_r8 .and. abs(maxval(cdf_diff)) < 1e-15_r8) then
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
beta = 2.0_r8

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 @@ -106,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 @@ -129,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 @@ -184,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 @@ -251,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 @@ -339,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 @@ -389,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