Created
May 12, 2024 17:42
-
-
Save stefanschmidt/d4cd1d8622c288a43b3af26ee82f01ee to your computer and use it in GitHub Desktop.
Compute the 75th percentile of the standard normal distribution with ndtri algorithm
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# Approximate the median of the standard half-normal distribution by | |
# computing the 75th percentile of the standard normal distribution using | |
# the algorithm from the original implementation of scipy.special.ndtri | |
# | |
# Based on the code from https://stackoverflow.com/questions/41338539#41338979 | |
# which was based on the original implementation in ndtri.c | |
s2pi = 2.50662827463100050242E0 | |
P0 = [ | |
-5.99633501014107895267E1, | |
9.80010754185999661536E1, | |
-5.66762857469070293439E1, | |
1.39312609387279679503E1, | |
-1.23916583867381258016E0, | |
] | |
Q0 = [ | |
1, | |
1.95448858338141759834E0, | |
4.67627912898881538453E0, | |
8.63602421390890590575E1, | |
-2.25462687854119370527E2, | |
2.00260212380060660359E2, | |
-8.20372256168333339912E1, | |
1.59056225126211695515E1, | |
-1.18331621121330003142E0, | |
] | |
def polevl(x, coef): | |
accum = 0 | |
for c in coef: | |
accum = x * accum + c | |
return accum | |
y0 = 0.75 | |
y = y0 | |
y = y - 0.5 | |
y2 = y * y | |
x = y + y * (y2 * polevl(y2, P0) / polevl(y2, Q0)) | |
x = x * s2pi | |
print(x) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment