Skip to content

Kolmogorov-Smirnov distribution #422

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

Merged
NAThompson merged 14 commits into from
Sep 4, 2020
Merged

Kolmogorov-Smirnov distribution #422

Changes from 1 commit
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
945b432
Kolmogorov-Smirnov distribution #421
evanmiller Aug 19, 2020
de347f1
Fix copyright notices [CI SKIP]
evanmiller Aug 19, 2020
59c4da1
Implement skewness for K-S distribution [CI SKIP]
evanmiller Aug 20, 2020
73f0dec
Implement kurtosis for the K-S distribution
evanmiller Aug 21, 2020
8ef4d3f
Fix "0.0" conversion errors with multiprecision types
evanmiller Aug 22, 2020
91fd5f6
Kolmogorov-Smirnov build fixes
evanmiller Aug 22, 2020
fad8917
Kolmogorov-Smirnov needs more BOOST_MATH_STD_USING
evanmiller Aug 22, 2020
ae6c2d8
Kolmogorov-Smirnov build fixes
evanmiller Aug 23, 2020
1599c1d
Require C++11 features for Kolmogorov-Smirnov tests
evanmiller Aug 23, 2020
dafb602
Fix double-to-float warnings with double literals
evanmiller Aug 24, 2020
07791a7
Kolmogorov-Smirnov user documentation [CI SKIP]
evanmiller Aug 30, 2020
c1565d1
Fix code comment (and kick off CI build)
evanmiller Aug 31, 2020
30e325a
Make suggested changes [CI SKIP]
evanmiller Sep 2, 2020
a1fa5e4
Tweak testing logic for K-S mode [CI SKIP]
evanmiller Sep 2, 2020
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
Prev Previous commit
Next Next commit
Fix double-to-float warnings with double literals 
Should fix Appveyor build
  • Loading branch information
@evanmiller
evanmiller committed Aug 31, 2020
commit dafb60258451a2f44415cba8fc4fbc002b192bff
28 changes: 12 additions & 16 deletions include/boost/math/distributions/kolmogorov_smirnov.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -120,16 +120,12 @@ namespace boost { namespace math {
namespace detail {
template <class RealType>
inline RealType kolmogorov_smirnov_quantile_guess(RealType p) {
RealType k;
// Choose a starting point for the Newton-Raphson iteration
if (p > 0.9) {
k = 1.8 - 5 * (1 - p);
} else if (p < 0.3) {
k = p + 0.45;
} else {
k = p + 0.3;
}
return k;
if (p > 0.9)
return RealType(1.8) - 5 * (1 - p);
if (p < 0.3)
return p + RealType(0.45);
return p + RealType(0.3);
Copy link
Collaborator

@NAThompson NAThompson Sep 2, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suppose it's low priority, but casting a double precision number like 0.3 to RealType is a little inelegant since it's not representable. Not really a big deal, but I'd use RealType(3)/RealType(10).

Copy link
Contributor Author

@evanmiller evanmiller Sep 2, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, the whole function is inelegant – it's basically me drawing three straight lines by hand to fit the curve for an initial guess. So I'm not too worried about epsilon-sized round-off errors.

}

// d/dk (theta2(0, 1/(2*k*k/M_PI))/sqrt(2*k*k*M_PI))
Expand All @@ -145,7 +141,7 @@ RealType kolmogorov_smirnov_pdf_small_x(RealType x, RealType n, const Policy&) {
return static_cast<RealType>(0);
}
while (1) {
delta = exp(-0.5*(i+0.5)*(i+0.5)*pi2/x2n) * ((i+0.5)*(i+0.5)*pi2 - x2n);
delta = exp(-RealType(i+0.5)*RealType(i+0.5)*pi2/(2*x2n)) * (RealType(i+0.5)*RealType(i+0.5)*pi2 - x2n);

if (delta == 0.0)
break;
Expand All @@ -169,7 +165,7 @@ inline RealType kolmogorov_smirnov_pdf_large_x(RealType x, RealType n, const Pol
RealType eps = policies::get_epsilon<RealType, Policy>();
int i = 1;
while (1) {
delta = 8*x*i*i*exp(-2.0*i*i*x*x*n);
delta = 8*x*i*i*exp(-2*i*i*x*x*n);

if (delta == 0.0)
break;
Expand All @@ -178,7 +174,7 @@ inline RealType kolmogorov_smirnov_pdf_large_x(RealType x, RealType n, const Pol
break;

if (i%2 == 0)
delta *= -1.0;
delta = -delta;

value += delta;
last_delta = delta;
Expand Down Expand Up @@ -447,8 +443,8 @@ inline RealType variance(const kolmogorov_smirnov_distribution<RealType, Policy>
RealType error_result;
if(false == detail::check_df(function, n, &error_result, Policy()))
return error_result;
return 0.5 * (constants::pi_sqr_div_six<RealType>()
- constants::pi<RealType>() * constants::ln_two<RealType>() * constants::ln_two<RealType>()) / n;
return (constants::pi_sqr_div_six<RealType>()
- constants::pi<RealType>() * constants::ln_two<RealType>() * constants::ln_two<RealType>()) / (2*n);
}

// Skewness and kurtosis come from integrating the PDF
Expand All @@ -462,7 +458,7 @@ inline RealType skewness(const kolmogorov_smirnov_distribution<RealType, Policy>
RealType error_result;
if(false == detail::check_df(function, n, &error_result, Policy()))
return error_result;
RealType ex3 = 0.5625 * constants::root_half_pi<RealType>() * constants::zeta_three<RealType>() / n / sqrt(n);
RealType ex3 = RealType(0.5625) * constants::root_half_pi<RealType>() * constants::zeta_three<RealType>() / n / sqrt(n);
RealType mean = boost::math::mean(dist);
RealType var = boost::math::variance(dist);
return (ex3 - 3 * mean * var - mean * mean * mean) / var / sqrt(var);
Expand All @@ -477,7 +473,7 @@ inline RealType kurtosis(const kolmogorov_smirnov_distribution<RealType, Policy>
RealType error_result;
if(false == detail::check_df(function, n, &error_result, Policy()))
return error_result;
RealType ex4 = 7.0 * constants::pi_sqr_div_six<RealType>() * constants::pi_sqr_div_six<RealType>() / 20.0 / n / n;
RealType ex4 = 7 * constants::pi_sqr_div_six<RealType>() * constants::pi_sqr_div_six<RealType>() / 20 / n / n;
RealType mean = boost::math::mean(dist);
RealType var = boost::math::variance(dist);
RealType skew = boost::math::skewness(dist);
Expand Down