Merge pull request #79 from sot/binom

Add function to compute binomial point percent function

chandra_aca/__init__.py

@@ -1,7 +1,7 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 from .transform import *
 
-__version__ = '4.26.1'
+__version__ = '4.27'
 
 
 def test(*args, **kwargs):

chandra_aca/star_probs.py

@@ -15,12 +15,9 @@
 
 """
 
-from __future__ import print_function, division
-
 import os
 import warnings
 from numba import jit
-from six.moves import zip
 
 from scipy.optimize import brentq, bisect
 import scipy.stats
@@ -726,3 +723,35 @@ def merit_func(t_ccd):
         return count - min_guide_count
 
     return bisect(merit_func, cold_t_ccd, warm_t_ccd, xtol=0.001, rtol=1e-15, full_output=False)
+
+
+def binom_ppf(k, n, conf, n_sample=1000):
+    """
+    Compute percent point function (inverse of CDF) for binomial, where
+    the percentage is with respect to the "p" (binomial probability) parameter
+    not the "k" parameter.  The latter is what one gets with scipy.stats.binom.ppf().
+
+    This function internally generates ``n_sample`` samples of the binomal PMF in
+    order to compute the CDF and finally interpolate to get the PPF.
+
+    The following example returns the 1-sigma (0.17 - 0.84) confidence interval
+    on the true binomial probability for an experiment with 4 successes in 5 trials.
+
+    Example::
+
+      >>> binom_ppf(4, 5, [0.17, 0.84])
+      array([ 0.55463945,  0.87748177])
+
+    :param k: int, number of successes (0 < k <= n)
+    :param n: int, number of trials
+    :param conf: float or array of floats, percent point values
+    :param n_sample: number of PMF samples for interpolation
+
+    :return: percent point function values corresponding to ``conf``
+    """
+    ps = np.linspace(0, 1, n_sample)  # Probability values
+    pmfs = scipy.stats.binom.pmf(k=k, n=n, p=ps)
+    cdf = np.cumsum(pmfs) / np.sum(pmfs)
+    out = np.interp(conf, xp=cdf, fp=ps)
+
+    return out

chandra_aca/tests/test_star_probs.py

@@ -9,7 +9,8 @@
 
 from chandra_aca.star_probs import (t_ccd_warm_limit, mag_for_p_acq, acq_success_prob,
                                     guide_count, t_ccd_warm_limit_for_guide,
-                                    grid_model_acq_prob, snr_mag_for_t_ccd)
+                                    grid_model_acq_prob, snr_mag_for_t_ccd,
+                                    binom_ppf)
 
 # Acquisition probabilities regression test data
 ACQ_PROBS_FILE = os.path.join(os.path.dirname(__file__), 'data', 'acq_probs.dat')
@@ -294,3 +295,8 @@ def test_grid_floor_2018_11():
                      -0.875, -0.695, -0.382, -0.204, 0.386, 0.758, 0.938, 1.133,
                      1.476, 1.639])
     assert np.allclose(probs.flatten(), exp, rtol=0, atol=0.001)
+
+
+def test_binom_ppf():
+    vals = binom_ppf(4, 5, [0.17, 0.84])
+    assert np.allclose(vals, [0.55463945, 0.87748177])