from math import erf, sqrt
from scipy.stats import norm # only for comparison
y = 0.123
z = 2*y-1
a = 0
while erf(a) > z or erf(a+1) < z: # looking for initial bracket of size 1
if erf(a) > z:
a -= 1
else:
a += 1
b = a+1 # found a bracket, proceed to refine it
while b-a > 1e-15: # 1e-15 ought to be enough precision
c = (a+b)/2.0 # bisection method
if erf(c) > z:
b = c
else:
a = c
print sqrt(2)*(a+b)/2.0 # this is the answer
print norm.ppf(y) # SciPy for comparison
没有一种简单的方法可以使用
erf
来实现norm.ppf
,因为norm.ppf
与erf
的逆有关。相反,这里是scipy
中代码的纯Python实现。您应该发现函数ndtri
返回的值与norm.ppf
完全相同:函数ppf是y=(1+erf(x/sqrt(2))/2的逆函数。所以我们需要解这个方程,给定y在0和1之间。下面是一个由bisection method执行此操作的代码。我导入了SciPy函数来说明结果是一样的。在
留给你去做:
相关问题 更多 >
编程相关推荐