Test that Bernoulli distribution has p=1 (i.e. non-random) from a sample

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$$?

• Welcome to CV. See our posts about the Rule of Three.
– whuber
Commented Oct 16, 2023 at 23:57
• @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$? Commented Oct 17, 2023 at 0:41
• Infinity is the answer. Commented Oct 17, 2023 at 10:55
• 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.
– whuber
Commented Oct 17, 2023 at 11:51
• He wants prove equality Commented Oct 17, 2023 at 14:12

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) } }