2
$\begingroup$

I have a black-box function $f$ that may only return 0 and 1 for a given input. The function may be deterministic or stochastic for different inputs.

Assume I have input $x_0$

How many times do I need to sample $f(x_0)$ to conclude that the function is deterministic with at least a probability $p_d$?

The problem is easy if a sample contains both $0$ and $1$ (this means that the function is non-determenistic), but what do I do if the sample is all $1$?

From my understanding of statistics, my function is essentially is a Bernoulli distribution sampler. How do I test for $p=1$ given a sample of all $1$?

$\endgroup$
6
  • 2
    $\begingroup$ Welcome to CV. See our posts about the Rule of Three. $\endgroup$
    – whuber
    Oct 16, 2023 at 23:57
  • $\begingroup$ @whuber I believe my problem is different. The rule of Three gives the lower bound for the $p$ with 95% confidence. I'm however interested in a little bit different question: how long do I need to sample $f(x_0)$ to make a conclusion that it is deterministic with a probability $p_d$? $\endgroup$ Oct 17, 2023 at 0:41
  • $\begingroup$ Infinity is the answer. $\endgroup$ Oct 17, 2023 at 10:55
  • $\begingroup$ The Rule of Three comes about because $-\log(1 - 95\%)\approx 3.$ It generalizes. @BigBendRegion appears to have a different interpretation, though -- but I don't know what it might be. $\endgroup$
    – whuber
    Oct 17, 2023 at 11:51
  • $\begingroup$ He wants prove equality $\endgroup$ Oct 17, 2023 at 14:12

1 Answer 1

1
$\begingroup$

To give this problem some notational structure, suppose we stipulate that we have an observable sequence $y_1,y_2,y_3, \cdots$ consisting of function outputs $f(x_0) \in \{ 0,1 \}$ using repeated inputs of $x_0$. We will use a Bernoulli model with $y_1,y_2,y_3, \cdots \sim \text{IID Bern}(\theta)$ with unknown probability $0 \leqslant \theta \leqslant 1$. Let $\mathscr{D}$ denote the hypothesis that the function is deterministic --- i.e.:

$$\mathscr{D} \quad \quad \iff \quad \quad y_1 = y_2 = y_3 = \cdots.$$

One possible approach here would be to use a Bayesian model with prior probability $\pi_0 = \mathbb{P}(\mathscr{D})$ and a beta prior $\theta | \bar{\mathscr{D}} \sim \text{Beta}(\alpha, \beta)$. (We will also assume that if the function is deterministic then there is an equal chance that it is deterministic to a zero or one.) Letting $s_n \equiv \sum_{i=1}^n y_i$ we then have:

$$\begin{align} p(y_1,...,y_n|\bar{\mathscr{D}}) &= \int \limits_0^1 p(y_1,...,y_n|\theta) \cdot p(\theta|\mathscr{D}) \ d\theta \\[6pt] &= \frac{1}{\text{B}(\alpha, \beta)} \int \limits_0^1 \theta^{\alpha+s_n-1}(1-\theta)^{\beta+n-s_n-1} \ d\theta \\[16pt] &= \frac{\text{B}(\alpha+s_n, \beta+n-s_n)}{\text{B}(\alpha, \beta)}, \\[22pt] p(y_1,...,y_n|\mathscr{D}) &= \frac{1}{2} \cdot \mathbb{I}(y_1 = \cdots = y_n). \\[6pt] \end{align}$$

where $\text{B}$ refers to the beta function. We then have the posterior probability:

$$\begin{align} \mathbb{P}(\mathscr{D} | y_1,...,y_n) &= \frac{p(y_1,...,y_n|\mathscr{D}) \cdot \mathbb{P}(\mathscr{D})}{p(y_1,...,y_n|\mathscr{D}) \cdot \mathbb{P}(\mathscr{D}) + p(y_1,...,y_n|\bar{\mathscr{D}}) \cdot \mathbb{P}(\bar{\mathscr{D}})} \\[12pt] &= \frac{\text{B}(\alpha, \beta) \cdot \pi_0}{\text{B}(\alpha, \beta) \cdot \pi_0 + 2 \text{B}(\alpha+s_n, \beta+n-s_n) \cdot (1-\pi_0)} \cdot \mathbb{I}(y_1 = \cdots = y_n), \\[12pt] \end{align}$$

which can be written in log-space (for $y_1 = \cdots = y_n$) as:

$$\begin{align} \log \mathbb{P}(\mathscr{D} | y_1,...,y_n) &= -\text{log1pexp} \bigg( \begin{matrix} \log \text{B}(\alpha+s_n, \beta+n-s_n) - \log \text{B}(\alpha, \beta) \\ + \log 2 + \log(1-\pi_0) - \log (\pi_0) \end{matrix} \bigg). \\[12pt] \end{align}$$

You can program the vectorised posterior probability function in R as follows:

POSTERIOR <- function(prior, n, s, alpha = 1, beta = 1, log = FALSE) {
  
  #Set log-beta functions
  LB1 <- lbeta(alpha, beta)
  LB2 <- lbeta(alpha+s, beta+n-s)
  
  #Compute posterior log-probability
  OUT <- rep(-Inf, length(prior))
  if ((s == 0)|(s == n)) {
    OUT <- -VGAM::log1pexp(LB2 - LB1 + log(2) + log(1-prior) - log(prior)) }
  
  #Return output
  if (log) { OUT } else { exp(OUT) } }
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.