def bayesianUpperLimit(self, alpha, steps=1.e5, plot=False):
"""
Compute one-sided upper limit using Bayesian Method of Helene.
Several methods of increasing numerical stability have been implemented.
"""
x_dense, y_dense = self.densify()
y_dense -= numpy.max(y_dense) # Numeric stability
f = scipy.interpolate.interp1d(x_dense, y_dense, kind='linear')
x = numpy.linspace(0., numpy.max(x_dense), steps)
pdf = numpy.exp(f(x) / 2.)
cut = (pdf / numpy.max(pdf)) > 1.e-10
x = x[cut]
pdf = pdf[cut]
#pdf /= pdf[0]
#forbidden = numpy.nonzero(pdf < 1.e-10)[0]
#if len(forbidden) > 0:
# index = forbidden[0] # Numeric stability
# x = x[0: index]
# pdf = pdf[0: index]
cdf = numpy.cumsum(pdf)
cdf /= cdf[-1]
cdf_reflect = scipy.interpolate.interp1d(cdf, x)
#if plot:
# pylab.figure()
# pylab.plot(x, f(x))
# pylab.scatter(self.x, self.y, c='red')
#
# pylab.figure()
# pylab.plot(x, pdf)
#
# pylab.figure()
# pylab.plot(cdf, x)
return cdf_reflect(alpha)
评论列表
文章目录