Skip to content

Instantly share code, notes, and snippets.

@tommyettinger
Created August 13, 2020 23:49
Show Gist options
  • Save tommyettinger/967b82f2e2f81f7a6929a3f9ce59abc1 to your computer and use it in GitHub Desktop.
Save tommyettinger/967b82f2e2f81f7a6929a3f9ce59abc1 to your computer and use it in GitHub Desktop.
Probit function in Java
/**
* Provides a {@link #probit(double)} method to convert from the (0.0,1.0) domain to a Gaussian range.
* Uses an algorithm by Peter John Acklam, as implemented by Sherali Karimov.
* Free to use, but you should credit Acklam and Karimov appropriately.
*/
public final class Probit {
private static final double LOW = 0.02425;
private static final double HIGH = 1.0 - LOW;
// Coefficients in rational approximations.
private static final double
A0 = -3.969683028665376e+01,
A1 = 2.209460984245205e+02,
A2 = -2.759285104469687e+02,
A3 = 1.383577518672690e+02,
A4 = -3.066479806614716e+01,
A5 = 2.506628277459239e+00;
private static final double
B0 = -5.447609879822406e+01,
B1 = 1.615858368580409e+02,
B2 = -1.556989798598866e+02,
B3 = 6.680131188771972e+01,
B4 = -1.328068155288572e+01;
private static final double
C0 = -7.784894002430293e-03,
C1 = -3.223964580411365e-01,
C2 = -2.400758277161838e+00,
C3 = -2.549732539343734e+00,
C4 = 4.374664141464968e+00,
C5 = 2.938163982698783e+00;
private static final double
D0 = 7.784695709041462e-03,
D1 = 3.224671290700398e-01,
D2 = 2.445134137142996e+00,
D3 = 3.754408661907416e+00;
/**
* A way of taking a double in the (0.0, 1.0) range and mapping it to a Gaussian or normal distribution, so high
* inputs correspond to high outputs, and similarly for the low range. This is centered on 0.0 and its standard
* deviation seems to be 1.0 (the same as {@link Random#nextGaussian()}). It is slightly faster (about a 10%
* increase in throughput, with comparable overhead for both methods) than the approach used by
* {@link Random#nextGaussian()}, even when both use a faster RNG algorithm and when both values produced by
* nextGaussian() are used.
* <br>
* This uses an algorithm by Peter John Acklam, as implemented by Sherali Karimov.
* <a href="https://web.archive.org/web/20150910002142/http://home.online.no/~pjacklam/notes/invnorm/impl/karimov/StatUtil.java">Source</a>.
* <a href="https://web.archive.org/web/20151030215612/http://home.online.no/~pjacklam/notes/invnorm/">Information on the algorithm</a>.
* @param d should be between 0 and 1, exclusive, but other values are tolerated (they return infinite results)
* @return a normal-distributed double centered on 0.0
*/
@SuppressWarnings("divzero") // This can legitimately return infinite doubles, which it produces with zero division.
public static double probit(final double d) {
if (d <= 0 || d >= 1) {
return (d - 0.5) / 0.0;
}
// Rational approximation for lower region:
else if (d < LOW) {
final double q = Math.sqrt(-2 * Math.log(d));
return (((((C0 * q + C1) * q + C2) * q + C3) * q + C4) * q + C5) / ((((D0 * q + D1) * q + D2) * q + D3) * q + 1);
}
// Rational approximation for upper region:
else if (HIGH < d) {
final double q = Math.sqrt(-2 * Math.log(1 - d));
return -(((((C0 * q + C1) * q + C2) * q + C3) * q + C4) * q + C5) / ((((D0 * q + D1) * q + D2) * q + D3) * q + 1);
}
// Rational approximation for central region:
else {
final double q = d - 0.5;
final double r = q * q;
return (((((A0 * r + A1) * r + A2) * r + A3) * r + A4) * r + A5) * q / (((((B0 * r + B1) * r + B2) * r + B3) * r + B4) * r + 1);
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment