def gpinv(p, k, sigma):
"""Inverse Generalised Pareto distribution function."""
x = np.empty(p.shape)
x.fill(np.nan)
if sigma <= 0:
return x
ok = (p > 0) & (p < 1)
if np.all(ok):
if np.abs(k) < np.finfo(float).eps:
np.negative(p, out=x)
np.log1p(x, out=x)
np.negative(x, out=x)
else:
np.negative(p, out=x)
np.log1p(x, out=x)
x *= -k
np.expm1(x, out=x)
x /= k
x *= sigma
else:
if np.abs(k) < np.finfo(float).eps:
# x[ok] = - np.log1p(-p[ok])
temp = p[ok]
np.negative(temp, out=temp)
np.log1p(temp, out=temp)
np.negative(temp, out=temp)
x[ok] = temp
else:
# x[ok] = np.expm1(-k * np.log1p(-p[ok])) / k
temp = p[ok]
np.negative(temp, out=temp)
np.log1p(temp, out=temp)
temp *= -k
np.expm1(temp, out=temp)
temp /= k
x[ok] = temp
x *= sigma
x[p == 0] = 0
if k >= 0:
x[p == 1] = np.inf
else:
x[p == 1] = -sigma / k
return x
评论列表
文章目录