1# Copyright (c) 2012, 2013 Ricardo Andrade 2# Copyright (c) 2015 James Hensman 3# Licensed under the BSD 3-clause license (see LICENSE.txt) 4 5import numpy as np 6from scipy.special import ndtr as std_norm_cdf 7 8#define a standard normal pdf 9_sqrt_2pi = np.sqrt(2*np.pi) 10def std_norm_pdf(x): 11 x = np.clip(x,-1e150,1e150) 12 return np.exp(-np.square(x)/2)/_sqrt_2pi 13 14def inv_std_norm_cdf(x): 15 """ 16 Inverse cumulative standard Gaussian distribution 17 Based on Winitzki, S. (2008) 18 """ 19 z = 2*x -1 20 ln1z2 = np.log(1-z**2) 21 a = 8*(np.pi -3)/(3*np.pi*(4-np.pi)) 22 b = 2/(np.pi * a) + ln1z2/2 23 inv_erf = np.sign(z) * np.sqrt( np.sqrt(b**2 - ln1z2/a) - b ) 24 return np.sqrt(2) * inv_erf 25 26def logPdfNormal(z): 27 """ 28 Robust implementations of log pdf of a standard normal. 29 30 @see [[https://github.com/mseeger/apbsint/blob/master/src/eptools/potentials/SpecfunServices.h original implementation]] 31 in C from Matthias Seeger. 32 """ 33 return -0.5 * (M_LN2PI + z * z) 34 35def cdfNormal(z): 36 """ 37 Robust implementations of cdf of a standard normal. 38 39 @see [[https://github.com/mseeger/apbsint/blob/master/src/eptools/potentials/SpecfunServices.h original implementation]] 40 in C from Matthias Seeger. 41 */ 42 """ 43 if (abs(z) < ERF_CODY_LIMIT1): 44 # Phi(z) approx (1+y R_3(y^2))/2, y=z/sqrt(2) 45 return 0.5 * (1.0 + (z / M_SQRT2) * _erfRationalHelperR3(0.5 * z * z)) 46 elif (z < 0.0): 47 # Phi(z) approx N(z)Q(-z)/(-z), z<0 48 return np.exp(logPdfNormal(z)) * _erfRationalHelper(-z) / (-z) 49 else: 50 return 1.0 - np.exp(logPdfNormal(z)) * _erfRationalHelper(z) / z 51 52 53 54def logCdfNormal(z): 55 """ 56 Robust implementations of log cdf of a standard normal. 57 58 @see [[https://github.com/mseeger/apbsint/blob/master/src/eptools/potentials/SpecfunServices.h original implementation]] 59 in C from Matthias Seeger. 60 """ 61 if (abs(z) < ERF_CODY_LIMIT1): 62 # Phi(z) approx (1+y R_3(y^2))/2, y=z/sqrt(2) 63 return np.log1p((z / M_SQRT2) * _erfRationalHelperR3(0.5 * z * z)) - M_LN2 64 elif (z < 0.0): 65 # Phi(z) approx N(z)Q(-z)/(-z), z<0 66 return logPdfNormal(z) - np.log(-z) + np.log(_erfRationalHelper(-z)) 67 else: 68 return np.log1p(-(np.exp(logPdfNormal(z))) * _erfRationalHelper(z) / z) 69 70 71 72def derivLogCdfNormal(z): 73 """ 74 Robust implementations of derivative of the log cdf of a standard normal. 75 76 @see [[https://github.com/mseeger/apbsint/blob/master/src/eptools/potentials/SpecfunServices.h original implementation]] 77 in C from Matthias Seeger. 78 """ 79 if (abs(z) < ERF_CODY_LIMIT1): 80 # Phi(z) approx (1 + y R_3(y^2))/2, y = z/sqrt(2) 81 return 2.0 * np.exp(logPdfNormal(z)) / (1.0 + (z / M_SQRT2) * _erfRationalHelperR3(0.5 * z * z)) 82 elif (z < 0.0): 83 # Phi(z) approx N(z) Q(-z)/(-z), z<0 84 return -z / _erfRationalHelper(-z) 85 else: 86 t = np.exp(logPdfNormal(z)) 87 return t / (1.0 - t * _erfRationalHelper(z) / z) 88 89 90def _erfRationalHelper(x): 91 assert x > 0.0, "Arg of erfRationalHelper should be >0.0; was {}".format(x) 92 93 if (x >= ERF_CODY_LIMIT2): 94 """ 95 x/sqrt(2) >= 4 96 97 Q(x) = 1 + sqrt(pi) y R_1(y), 98 R_1(y) = poly(p_j,y) / poly(q_j,y), where y = 2/(x*x) 99 100 Ordering of arrays: 4,3,2,1,0,5 (only for numerator p_j; q_5=1) 101 ATTENTION: The p_j are negative of the entries here 102 p (see P1_ERF) 103 q (see Q1_ERF) 104 """ 105 y = 2.0 / (x * x) 106 107 res = y * P1_ERF[5] 108 den = y 109 i = 0 110 111 while (i <= 3): 112 res = (res + P1_ERF[i]) * y 113 den = (den + Q1_ERF[i]) * y 114 i += 1 115 116 # Minus, because p(j) values have to be negated 117 return 1.0 - M_SQRTPI * y * (res + P1_ERF[4]) / (den + Q1_ERF[4]) 118 else: 119 """ 120 x/sqrt(2) < 4, x/sqrt(2) >= 0.469 121 122 Q(x) = sqrt(pi) y R_2(y), 123 R_2(y) = poly(p_j,y) / poly(q_j,y), y = x/sqrt(2) 124 125 Ordering of arrays: 7,6,5,4,3,2,1,0,8 (only p_8; q_8=1) 126 p (see P2_ERF) 127 q (see Q2_ERF 128 """ 129 y = x / M_SQRT2 130 res = y * P2_ERF[8] 131 den = y 132 i = 0 133 134 while (i <= 6): 135 res = (res + P2_ERF[i]) * y 136 den = (den + Q2_ERF[i]) * y 137 i += 1 138 139 return M_SQRTPI * y * (res + P2_ERF[7]) / (den + Q2_ERF[7]) 140 141def _erfRationalHelperR3(y): 142 assert y >= 0.0, "Arg of erfRationalHelperR3 should be >=0.0; was {}".format(y) 143 144 nom = y * P3_ERF[4] 145 den = y 146 i = 0 147 while (i <= 2): 148 nom = (nom + P3_ERF[i]) * y 149 den = (den + Q3_ERF[i]) * y 150 i += 1 151 return (nom + P3_ERF[3]) / (den + Q3_ERF[3]) 152 153ERF_CODY_LIMIT1 = 0.6629 154ERF_CODY_LIMIT2 = 5.6569 155M_LN2PI = 1.83787706640934533908193770913 156M_LN2 = 0.69314718055994530941723212146 157M_SQRTPI = 1.77245385090551602729816748334 158M_SQRT2 = 1.41421356237309504880168872421 159 160#weights for the erfHelpers (defined here to avoid redefinitions at every call) 161P1_ERF = [ 1623.05326634961232344e-1, 3.60344899949804439e-1, 1631.25781726111229246e-1, 1.60837851487422766e-2, 1646.58749161529837803e-4, 1.63153871373020978e-2] 165Q1_ERF = [ 1662.56852019228982242e+0, 1.87295284992346047e+0, 1675.27905102951428412e-1, 6.05183413124413191e-2, 1682.33520497626869185e-3] 169P2_ERF = [ 1705.64188496988670089e-1, 8.88314979438837594e+0, 1716.61191906371416295e+1, 2.98635138197400131e+2, 1728.81952221241769090e+2, 1.71204761263407058e+3, 1732.05107837782607147e+3, 1.23033935479799725e+3, 1742.15311535474403846e-8] 175Q2_ERF = [ 1761.57449261107098347e+1, 1.17693950891312499e+2, 1775.37181101862009858e+2, 1.62138957456669019e+3, 1783.29079923573345963e+3, 4.36261909014324716e+3, 1793.43936767414372164e+3, 1.23033935480374942e+3] 180P3_ERF = [ 1813.16112374387056560e+0, 1.13864154151050156e+2, 1823.77485237685302021e+2, 3.20937758913846947e+3, 1831.85777706184603153e-1] 184Q3_ERF = [ 1852.36012909523441209e+1, 2.44024637934444173e+2, 1861.28261652607737228e+3, 2.84423683343917062e+3] 187