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.

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

bug: beta distribution not correct (failing tests) #717

Open
hkershaw-brown opened this issue Aug 13, 2024 · 9 comments
Open

bug: beta distribution not correct (failing tests) #717

hkershaw-brown opened this issue Aug 13, 2024 · 9 comments
Labels
Bug Something isn't working QCEFF quantile conserving filters

Comments

@hkershaw-brown
Copy link
Member

🐛 Your bug may already be reported!
Please search on the issue tracker before creating a new issue.

Describe the bug

  1. Run https://github.com/NCAR/DART/tree/stats-tests/developer_tests/beta_dist
  2. What was the expected outcome? test passes
  3. 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:

  1. Mac M2
  2. GNU Fortran (MacPorts gcc13 13.2.0_4+stdlib_flag) 13.2.0
@hkershaw-brown hkershaw-brown added the Bug Something isn't working label Aug 13, 2024
@hkershaw-brown hkershaw-brown changed the title bug: beta distribution failing tests bug: beta distribution not correct (failing tests) Aug 14, 2024
@hkershaw-brown
Copy link
Member Author

5b3850f

@hkershaw-brown
Copy link
Member Author

https://github.com/NCAR/DART/tree/beta_distribution_fix

note: default values for distribution type NO_DISTRIBUTION, not bounded?

type distribution_params_type
integer :: distribution_type
logical :: bounded_below, bounded_above
real(r8) :: lower_bound, upper_bound
real(r8) :: params(2)
integer :: ens_size
real(r8), allocatable :: ens(:)
real(r8), allocatable :: more_params(:)
end type

@jlaucar
Copy link
Contributor

jlaucar commented Sep 12, 2024 via email

@hkershaw-brown
Copy link
Member Author

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.

@jlaucar
Copy link
Contributor

jlaucar commented Sep 12, 2024 via email

@hkershaw-brown
Copy link
Member Author

hkershaw-brown commented Sep 13, 2024

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

ba3513e

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

@hkershaw-brown
Copy link
Member Author

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

@hkershaw-brown
Copy link
Member Author

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
 --------------------------------------

@hkershaw-brown
Copy link
Member Author

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Bug Something isn't working QCEFF quantile conserving filters
Projects
None yet
Development

No branches or pull requests

2 participants