Created
March 15, 2022 15:41
-
-
Save zimolzak/dec914b336c53f7c93b31ee3f6327d59 to your computer and use it in GitHub Desktop.
Inverse survival function of standard normal: pretty-print the polynomial approximation
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
print("""What is the "inverse survival function" of the standard normal? | |
The cdf is based on erf(x), the survival function is based on erfc(x), | |
and the inverse survival is based on 1/erfc(x). The error function | |
erf(x) is not an elementary function, but it can be approximated as | |
1 - 1 / (p(x))^4 , where p(x) is a 4th degree polynomial in x, with | |
coefficients a = 0.278393, b = 0.230389, c = 0.000972, d = 0.078108. | |
[https://en.wikipedia.org/wiki/Error_function] This means that | |
1/erfc(x) can be approximated as (p(x))^4 . | |
Specifically, (1 + a x + b x^2 + c x^3 + d x^4)^4 . | |
Wolfram|Alpha will expand that for us, but it's not nicely organized | |
by exponent of x. | |
Hence : | |
""") | |
wolfram_output = "a^4 x^4 + 4 a^3 b x^5 + 4 a^3 c x^6 + 4 a^3 d x^7 + 4 a^3 x^3 + 6 a^2 b^2 x^6 + 12 a^2 b c x^7 + 12 a^2 b d x^8 + 12 a^2 b x^4 + 6 a^2 c^2 x^8 + 12 a^2 c d x^9 + 12 a^2 c x^5 + 6 a^2 d^2 x^10 + 12 a^2 d x^6 + 6 a^2 x^2 + 4 a b^3 x^7 + 12 a b^2 c x^8 + 12 a b^2 d x^9 + 12 a b^2 x^5 + 12 a b c^2 x^9 + 24 a b c d x^10 + 24 a b c x^6 + 12 a b d^2 x^11 + 24 a b d x^7 + 12 a b x^3 + 4 a c^3 x^10 + 12 a c^2 d x^11 + 12 a c^2 x^7 + 12 a c d^2 x^12 + 24 a c d x^8 + 12 a c x^4 + 4 a d^3 x^13 + 12 a d^2 x^9 + 12 a d x^5 + 4 a x + b^4 x^8 + 4 b^3 c x^9 + 4 b^3 d x^10 + 4 b^3 x^6 + 6 b^2 c^2 x^10 + 12 b^2 c d x^11 + 12 b^2 c x^7 + 6 b^2 d^2 x^12 + 12 b^2 d x^8 + 6 b^2 x^4 + 4 b c^3 x^11 + 12 b c^2 d x^12 + 12 b c^2 x^8 + 12 b c d^2 x^13 + 24 b c d x^9 + 12 b c x^5 + 4 b d^3 x^14 + 12 b d^2 x^10 + 12 b d x^6 + 4 b x^2 + c^4 x^12 + 4 c^3 d x^13 + 4 c^3 x^9 + 6 c^2 d^2 x^14 + 12 c^2 d x^10 + 6 c^2 x^6 + 4 c d^3 x^15 + 12 c d^2 x^11 + 12 c d x^7 + 4 c x^3 + d^4 x^16 + 4 d^3 x^12 + 6 d^2 x^8 + 4 d x^4 + 1" | |
terms = wolfram_output.split(" + ") | |
MAX_EXPONENT = 16 | |
# Initialize empty list of lists. | |
polynomial = [] | |
for i in range(MAX_EXPONENT + 1): | |
polynomial += [[]] # [[]] * 17 won't work right | |
# Populate the list of lists, carefully keeping things ordered by the | |
# degree of each term. | |
for t in terms: | |
if not 'x' in t: | |
coefficient_str = t | |
degree_str = '0' | |
elif not '^' in t: | |
coefficient_str, degree_str = t.split(' x') | |
degree_str = '1' # reassign | |
else: | |
coefficient_str, degree_str = t.split(' x^') | |
degree = int(degree_str) | |
polynomial[degree] += [coefficient_str] | |
# Print the polynomial, in order of descending powers of x. | |
for degree in range(MAX_EXPONENT, -1, -1): | |
string_list = polynomial[degree] | |
coefficient_str = ' + '.join(string_list) | |
if degree == 0: | |
print(coefficient_str) # Just say the constant, not x^0. | |
elif degree == 1: | |
print("(%s) x +" % (coefficient_str)) # Just say x, not x^1. | |
else: | |
print("(%s) x^%i +" % (coefficient_str, degree)) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment