NCAR
bug: beta distribution not correct (failing tests)
:bug: Your bug may already be reported! Please search on the issue tracker before creating a new issue.
Describe the bug
- Run https://github.com/NCAR/DART/tree/stats-tests/developer_tests/beta_dist
- What was the expected outcome? test passes
- What actually happened? test fails: max difference in inversion is 7.8846335608728112E-006 max difference should be less than 1e-14
Error Message
[hkershaw:work] (stats-tests) > ./test_beta_dist
--------------------------------------
Starting ... at YYYY MM DD HH MM SS =
2024 8 13 15 13 6
--------------------------------------
set_nml_output Echo NML values to log file only
Absolute value of differences should be less than 1e-15
1 -2.2204460492503131E-016 3.0531133177191805E-016
2 -6.9388939039072284E-018 -3.0357660829594124E-018
3 0.0000000000000000 0.0000000000000000
4 0.0000000000000000 -1.1102230246251565E-016
5 8.3266726846886741E-017 0.0000000000000000
6 0.0000000000000000 -1.1102230246251565E-016
7 -1.1102230246251565E-016 -2.2204460492503131E-016
inv_cdf Failed to converge for quantile 0.0000000000000000
----------------------------
max difference in inversion is 7.8846335608728112E-006
max difference should be less than 1e-14
--------------------------------------
Finished ... at YYYY MM DD HH MM SS =
2024 8 13 15 13 6
--------------------------------------
Which model(s) are you working with?
none
Version of DART
Which version of DART are you using? v11.6.0-4-g2feb86983
Have you modified the DART code?
Yes, added test harnesses https://github.com/NCAR/DART/tree/stats-tests
Build information
Please describe:
- Mac M2
- GNU Fortran (MacPorts gcc13 13.2.0_4+stdlib_flag) 13.2.0
https://github.com/NCAR/DART/commit/5b3850fb1a08fff0ebce3f95283fef29de431bdc
https://github.com/NCAR/DART/tree/beta_distribution_fix
note: default values for distribution type NO_DISTRIBUTION, not bounded?
https://github.com/NCAR/DART/blob/71d70e1e346e3cc7c0244ad839415104140cfc9d/assimilation_code/modules/assimilation/distribution_params_mod.f90#L10-L18
The normal distribution should be unbounded.
On Thu, Sep 12, 2024 at 1:42 PM Helen Kershaw @.***> wrote:
https://github.com/NCAR/DART/tree/beta_distribution_fix
note: default values for distribution type NO_DISTRIBUTION, not bounded?
https://github.com/NCAR/DART/blob/71d70e1e346e3cc7c0244ad839415104140cfc9d/assimilation_code/modules/assimilation/distribution_params_mod.f90#L10-L18
— Reply to this email directly, view it on GitHub https://github.com/NCAR/DART/issues/717#issuecomment-2347102799, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANDHUIVEHBGTTKJ64NYALWTZWHVDDAVCNFSM6AAAAABMO4S2W2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGNBXGEYDENZZHE . You are receiving this because you are subscribed to this thread.Message ID: @.***>
Yeah I was just wondering if distribution_params_type if it should be initialized to distribution_type=NOT_A_DISTRIBUTION
so if it did not get set, it would trap on anything checking for known distributions.
Helen, Yes, doing that would probably be a safe practice. I have not tried to implement that or even comprehensively check for other distributions at this point. Jeff
On Thu, Sep 12, 2024 at 2:03 PM Helen Kershaw @.***> wrote:
Yeah I was just wondering if distribution_params_type if it should be initialized to distribution_type=NOT_A_DISTRIBUTION
so if it did not get set, it would trap on anything checking for known distributions.
— Reply to this email directly, view it on GitHub https://github.com/NCAR/DART/issues/717#issuecomment-2347134539, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANDHUITRDF4XPYHFGCFTVJDZWHXPJAVCNFSM6AAAAABMO4S2W2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGNBXGEZTINJTHE . You are receiving this because you commented.Message ID: @.***>
nest the if statement fix on https://github.com/NCAR/DART/tree/stats-fixes
This fix only works for gfortran13, not for ifort (IFORT) 2021.10.0
fix: split if statement, the compound logical condition does not get entered and gives a 'Failed to converge for quantile' when quantile==0
https://github.com/NCAR/DART/commit/ba3513e47ab2d2ba63c50137c897323ba0689836
Q. Is this a floating point comparison error? (plus short circuting the if statement)
-Wcompare-reals (for distribution mods only)
[hkershaw:work](main) > make
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/utilities/types_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/utilities/sort_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/utilities/utilities_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/distribution_params_mod.f90
mpif90 -Wcompare-reals -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/normal_distribution_mod.f90
/Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/normal_distribution_mod.f90:363:22:
363 | if(bounded_below .and. quantile == 0.0_r8) then
| 1
Warning: Equality comparison for REAL(8) at (1) [-Wcompare-reals]
/Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/normal_distribution_mod.f90:367:22:
367 | if(bounded_above .and. quantile == 1.0_r8) then
| 1
Warning: Equality comparison for REAL(8) at (1) [-Wcompare-reals]
/Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/normal_distribution_mod.f90:422:9:
422 | if(q_old == q_new) exit
| 1
Warning: Equality comparison for REAL(8) at (1) [-Wcompare-reals]
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/utilities/random_seq_mod.f90
mpif90 -Wcompare-reals -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/gamma_distribution_mod.f90
/Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/gamma_distribution_mod.f90:242:3:
242 | if(x == 0.0_r8) then
| 1
Warning: Equality comparison for REAL(8) at (1) [-Wcompare-reals]
/Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/gamma_distribution_mod.f90:294:10:
294 | if (pn(6) /= 0.0_r8) then
| 1
Warning: Inequality comparison for REAL(8) at (1) [-Wcompare-reals]
/Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/gamma_distribution_mod.f90:185:7:
185 | elseif(x == 0.0_r8) then
| 1
Warning: Equality comparison for REAL(8) at (1) [-Wcompare-reals]
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/beta_distribution_mod.f90
mpif90 -Wcompare-reals -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/bnrh_distribution_mod.f90
/Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/bnrh_distribution_mod.f90:621:9:
621 | if(sorted_ens(i) == lower_bound) then
| 1
Warning: Equality comparison for REAL(8) at (1) [-Wcompare-reals]
/Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/bnrh_distribution_mod.f90:633:9:
633 | if(sorted_ens(i) == upper_bound) then
| 1
Warning: Equality comparison for REAL(8) at (1) [-Wcompare-reals]
/Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/bnrh_distribution_mod.f90:652:6:
652 | if(sorted_ens(i) == sorted_ens(i - 1)) then
| 1
Warning: Equality comparison for REAL(8) at (1) [-Wcompare-reals]
/Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/bnrh_distribution_mod.f90:289:7:
289 | elseif(x == sort_ens(1)) then
| 1
Warning: Equality comparison for REAL(8) at (1) [-Wcompare-reals]
/Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/bnrh_distribution_mod.f90:335:13:
335 | elseif(x == sort_ens(j+1)) then
| 1
Warning: Equality comparison for REAL(8) at (1) [-Wcompare-reals]
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/utilities/time_manager_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/utilities/mpi_utilities_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/rootfinding_mod.f90
mpif90 -Wcompare-reals -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/kde_distribution_mod.f90
/Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/kde_distribution_mod.f90:121:16:
121 | if (x == 0._r8) then
| 1
Warning: Equality comparison for REAL(8) at (1) [-Wcompare-reals]
/Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/kde_distribution_mod.f90:144:16:
144 | if (x .eq. 0._r8) then
| 1
Warning: Equality comparison for REAL(8) at (1) [-Wcompare-reals]
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/probit_transform_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/utilities/ensemble_manager_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/location/utilities//default_location_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/location/oned/location_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/observations/obs_kind_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/utilities/netcdf_utilities_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/location/utilities//location_io_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/utilities/no_cray_win_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/utilities/distributed_state_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/io/state_structure_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/io/dart_time_io_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/models/utilities//default_model_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/models/lorenz_96/model_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/assim_model_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/cov_cutoff_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main//observations/forward_operators/obs_def_utilities_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/observations/forward_operators/obs_def_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/observations/obs_sequence_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/quality_control_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/observations/forward_operator_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/adaptive_inflate_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/utilities/options_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/io/io_filenames_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/io/direct_netcdf_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/io/state_vector_io_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/reg_factor_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/models/model_mod_tools/model_check_utilities_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/models/model_mod_tools/test_interpolate_oned.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/utilities/assert_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/models/utilities//quad_utils_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/algorithm_info_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/utilities/parse_args_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/sampling_error_correction_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/utilities/schedule_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/utilities/obs_impact_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/assim_tools_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/obs_model_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/modules/assimilation/filter_mod.f90
mpif90 -O2 -ffree-line-length-none -I/opt/local/include -c /Users/hkershaw/DART/Projects/NumericalComp/DART.main/assimilation_code/programs/filter/filter.f90
mpif90 probit_transform_mod.o types_mod.o forward_operator_mod.o ensemble_manager_mod.o gamma_distribution_mod.o dart_time_io_mod.o netcdf_utilities_mod.o bnrh_distribution_mod.o random_seq_mod.o rootfinding_mod.o state_vector_io_mod.o obs_kind_mod.o reg_factor_mod.o direct_netcdf_mod.o no_cray_win_mod.o quality_control_mod.o test_interpolate_oned.o state_structure_mod.o assert_mod.o cov_cutoff_mod.o quad_utils_mod.o obs_def_utilities_mod.o options_mod.o algorithm_info_mod.o parse_args_mod.o sampling_error_correction_mod.o obs_def_mod.o schedule_mod.o model_check_utilities_mod.o distribution_params_mod.o adaptive_inflate_mod.o normal_distribution_mod.o sort_mod.o io_filenames_mod.o assim_model_mod.o filter.o default_location_mod.o beta_distribution_mod.o model_mod.o distributed_state_mod.o mpi_utilities_mod.o location_io_mod.o time_manager_mod.o obs_model_mod.o utilities_mod.o assim_tools_mod.o default_model_mod.o location_mod.o obs_impact_mod.o kde_distribution_mod.o obs_sequence_mod.o filter_mod.o -o filter -O2 -ffree-line-length-none -I/opt/local/include -L/opt/local/lib -lnetcdff -lnetcdf
At least some of the "Failed to converge" seems to be the root finding oscillating in inv_cdf
quantile = 0.99999996152378834_dp
del_q = -1.1102230246251565e-16_dp
x_guess = 5.3741321112704998_dp
q_guess = 0.99999996152378934_dp
flips +/- del_q, del_q old over and over.
del_q = -1.1102230246251565e-16_dp
del_q_old = 1.1102230246251565e-16_dp
del_q = 1.1102230246251565e-16_dp
del_q_old = -1.1102230246251565e-16_dp
reproducer:
hkershaw@derecho4:/glade/derecho/scratch/hkershaw/DART/Bugs/beta_dist/DART/developer_tests/normal_dist/work(stats-fixes)$ cat ../test_oscillation.f90
program test_oscillation
use normal_distribution_mod, only : inv_std_normal_cdf
use types_mod, only : r4, r8, digits12
use utilities_mod, only: initialize_utilities, finalize_utilities
implicit none
real(r8) :: quantile
real(r8) :: x
call initialize_utilities()
quantile = 0.99999996152378834_r8
x = inv_std_normal_cdf(quantile)
call finalize_utilities()
end program test_oscillation
printout in normal_distribution_mod.f90
diff --git a/assimilation_code/modules/assimilation/normal_distribution_mod.f90 b/assimilation_code/modules/assimilation/normal_distribution_mod.f90
index f915ad66b..147bf6e62 100644
--- a/assimilation_code/modules/assimilation/normal_distribution_mod.f90
+++ b/assimilation_code/modules/assimilation/normal_distribution_mod.f90
@@ -418,9 +418,11 @@ do iter = 1, max_iterations
! If we've gone too far, the new error will be bigger than the old;
! Repeatedly half the distance until this is rectified
del_q_old = del_q
+ print*, 'del_q_old', del_q_old
q_new = cdf(x_new, p)
do j = 1, max_half_iterations
del_q = q_new - quantile
+ print*, 'del_q ', del_q
if (abs(del_q) < abs(del_q_old)) then
exit
endif
q_old = q_new
x_new = (x_guess + x_new)/2.0_r8
q_new = cdf(x_new, p)
! If q isn't changing, no point in continuing
if(q_old == q_new) exit
end do
output:
--------------------------------------
Starting ... at YYYY MM DD HH MM SS =
2024 10 4 13 4 9
--------------------------------------
set_nml_output Echo NML values to log file only
del_q_old 9.992007221626409E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
del_q_old -1.110223024625157E-016
del_q 1.110223024625157E-016
del_q_old 1.110223024625157E-016
del_q -1.110223024625157E-016
inv_cdf Failed to converge for quantile 0.999999961523788
--------------------------------------
Finished ... at YYYY MM DD HH MM SS =
2024 10 4 13 4 9
--------------------------------------
looking at the guesses for an oscillating case:
.....
x_new 5.37413211153138
x_guess 5.37413211153138
x_new 5.37413211127050
inv_cdf Failed to converge for quantile 0.999999961523788
x 5.37413211127050
If you do the normcdf in Matlab for both of these you get the same answer. If you do the norminv, you don't get either of them back.
>> normcdf(5.37413211153138)
ans =
0.999999961523788
>> normcdf(5.37413211127050)
ans =
0.999999961523788
>> norminv(0.999999961523788)
ans =
5.374132109841061
>> normcdf(5.374132109841061)
ans =
0.999999961523788
tests are in release v11.8.6
Fix is not released https://github.com/NCAR/DART/tree/beta_distribution_fix notes in this commit message: https://github.com/NCAR/DART/commit/5b3850fb1a08fff0ebce3f95283fef29de431bdc
For completeness, commit message from https://github.com/NCAR/DART/tree/beta_distribution_fix notes in this commit message: https://github.com/NCAR/DART/commit/5b3850fb1a08fff0ebce3f95283fef29de431bdc (branch will get deleted)
Commit 5b3850f jlaucar on Aug 20 The commit with sha cbe0ca on 5 April 2023 transitioned to the use of 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. beta_distribution_fix
Gamma_distribution_mod also changed in https://github.com/NCAR/DART/commit/5b3850fb1a08fff0ebce3f95283fef29de431bdc to only support supports the standard gamma distribution that is bounded below by 0 and unbounded above.
Q. what does this mean for user options (bounds) set in QCEFF table? Ignored for Gamma and Beta?