in bayesmark/quantiles.py [0:0]
def max_quantile_CI(X, q, m, alpha=0.05):
"""Calculate CI on `q` quantile of distribution on max of `m` iid samples using a data set `X`.
This uses nonparametric estimation from order statistics and will have alpha level of at most `alpha` due to the
discrete nature of order statistics.
Parameters
----------
X : :class:`numpy:numpy.ndarray` of shape (n,)
Data for quantile estimation. Can be vectorized. Must be sortable data type (which is almost everything).
q : float
Quantile to compute, must be in (0, 1). Can be vectorized.
m : int
Compute statistics for distribution on max over `m` samples. Must be ``>= 1``. Can be vectorized.
alpha : float
False positive rate we allow for CI, must be in (0, 1). Can be vectorized.
Returns
-------
estimate : dtype of `X`, scalar
Best estimate on `q` quantile on max over `m` iid samples.
LB : dtype of `X`, scalar
Lower end on CI
UB : dtype of `X`, scalar
Upper end on CI
"""
# X and alpha used/checked below in quantile_CI routine.
# We could robustify things to allow the edge cases, but maybe later
assert np.all(0 < q) and np.all(q < 1)
# Could check int but if someone wants to interpolate, we will let them.
assert np.all(m >= 1)
# Currently don't support broadcasting both at same time
assert np.ndim(X) == 1 or (np.ndim(q) == 0 and np.ndim(q) == 0 and np.ndim(alpha) == 0)
q = q ** (1.0 / m)
o_stats = order_stats(X)
n = X.shape[-1]
idx = _quantile(n, q)
idx_lower, idx_upper = _quantile_CI(n, q, alpha=alpha)
LB, UB = o_stats[..., idx_lower], o_stats[..., idx_upper]
# Might need to broadcast estimate out if vectorization is in alpha
estimate = ensure_shape(o_stats[..., idx], LB)
return estimate, LB, UB