Confidence intervals for a polynomial

Question

I have a random variable $Z$ which takes values in the nonnegative integers $\{ 0,1,2,\dots \}$, call the probabilities for each outcome $z_k:=P[Z=k]$. I can sample from $Z$'s distribution independently and cheaply; I currently have a sample size of $2^{28}$. It appears that $z_0\approx 0.24, z_1\approx 0.18,\dots$, with roughly exponential decay.

I have a sequence of quadratic forms with positive coefficients:

• $Q_0(z_0) = \frac14 z_0^2$
• $Q_1(z_0,z_1) = \frac 12 {z_0 z_1}$
• ...
• $Q_7(z_0,z_1,\dots,z_7) = \frac{1}{8} \left(2 z_0 z_1+3 z_2 z_1+4 z_4 z_1+4 z_6 z_1+3 z_0 z_3 + \right.$ $\qquad\qquad\qquad\qquad\qquad\qquad \left. +4 z_2 z_3+4 z_3 z_4+4 z_0 z_5+4 z_2 z_5+4 z_0 z_7\right)$
• ...

What I would like to have is a confidence interval for the $Q_i$'s that is less than $10^{-4}$ wide, but I'll take whatever I can get.

I have rigorous bounds on the $z_i$, and since the coefficients of the $Q$'s are all positive, it is straightforward to turn these into rigorous bounds for the $Q$'s. But I don't know how to do this correctly with confidence intervals.

---

What's this about? I found a bizarre phenomenon in number theory, and I know how to prove that it really happens, but actually doing so will require some programming effort on my part and a considerable amount of time on our local cluster. Before I invest that time and clog up our machine, I'd like to be more certain than I am that the phenomenon is real.

I want to quantify the reasonableness of my claim that $Q_7<Q_6$ and $Q_7<Q_8$. My estimates indicate that $Q_6-Q_7$ is around $5\cdot 10^{-4}$, which is why I wanted CIs at that resolution.

Fix a large integer $n$, and let $A$ be a uniformly chosen subset of $\{1,2,\dots,n\}$ (that is, each particular subset has probability $2^{-n}$ of being chosen). Let $Q_k(n)$ be the probability that exactly $k$ of the numbers from $\{2,3,\dots,2n\}$ cannot be written as a sum of two elements of $A$; let $Q_k = \lim_n Q_k(n)$. It's a little tricky to prove, but those limits exist and $\sum_{k} Q_k =1$. Now it's no surprise that $Q_0$ is small, and as $k$ increases $Q_k$ increases, has a peak and then decays exponentially. The bizarre part is that there is a bias against 7. That is, experimentally $Q_7< Q_6$ and $Q_7<Q_8$. That is, what wasn't a surprise actually isn't true: the distribution is bimodal.

I can express the $Q_i$'s (using some theory) as above without the limit in terms of this other distribution, defined by the $z_i$'s. That's handy because I have a way to rigorously bound the $z_i$'s using, as I mentioned above, some large computations. Also, I have a very large data set for the $Z$ variable.

Answer 1

In my answer, I provide many links to background material to save space here. I'm going to write my answer taking the info in the links as given.

I think a Bayesian approach is a natural fit to this problem, especially since you seek to convince only yourself. It's a bit convoluted to use confidence intervals to answer the question you really care about, to wit, how plausible is it that $Q_{7}<Q_{6}$ and $Q_{7}<Q_{8}$ given the sample from the $z_{i}$ distribution? The Bayesian approach allows you to address this question directly.

Likelihood function

Let $f_k$ be the observed frequency of integer outcome $k$ in your sample and let $N$ be the sample size. The likelihood function is proportional to the multinomial distribution. It has the form

$L(z_{0},...z_{8};f_{0},...f_{8})=\prod_{i=0}^{8}{z_{i}}^{Nf_{i}}$.

Prior distribution

The Dirichlet distribution is the natural choice for the prior distribution because it is the conjugate prior for the multinomial likelihood. It has the form

$p(z_{0},...z_{8};\alpha_{0},...,\alpha_{8})\propto\prod_{i=0}^{8}{z_{i}}^{\alpha_{i}-1}$

This prior has nine hyperparameters (the $\alpha_i$ values), and they're a bit of a pain to deal with. In this “large sample” context, any reasonable choice of hyperparameter values will have negligible influence on the result, but still, I think it's worth devoting a bit of effort to selecting sensible values.

Here's how I recommend setting the hyperparameters. First, note that under this distribution $\mathrm{E}(z_{i})=\frac{\alpha_{i}}{\sum_{i=0}^{8}\alpha_{i}}$. Next, note that the simplest maximum entropy distribution over the naturals is the geometric distribution. So set

$\alpha_{i+1}=r\alpha_{i}=r^{i}\alpha_{0},\,0<r<1,$

$\alpha_{0}=A\left(\frac{1-r}{1-r^{9}}\right).$

Then $\mathrm{E}(z_{i})=r^{i}\left(\frac{1-r}{1-r^{9}}\right)$, so the distribution of the $z_{i}$ values is centered on a (truncated) geometric distribution. Furthermore, $\mathrm{Var}\left(z_{i}\right)\propto\frac{1}{(A+1)}$, so the value of $A$ controls the dispersion around this expectation but has no effect on the expectation itself.

This specification reduces the number of hyperparameters from the nine $\alpha_{i}$ values to just $r$ and $A$. I'll defer discussion of specific values of $r$ and $A$ for now.

Posterior probability of the proposition of interest

The posterior distribution of the $z_{i}$ values is the following Dirichlet distribution:

$p(z_{0},...z_{8}|f_{0},...,f_{8})\propto\prod_{i=0}^{8}{z_{i}}^{\alpha_{i}+Nf_{i}-1}.$

Let $\mathbb{Y}=\left\{ z_{0},...z_{8}|Q_7<Q_6 \text{ and } Q_7<Q_8\right\} $. The posterior probability you're interested in is

$\Pr(Q_7<Q_6 \text{ and } Q_7<Q_8|f_0,...,f_8) \propto \int_{\mathbb{Y}}\prod_{i=0}^{8}{z_{i}}^{\alpha_{i}+Nf_i-1}dz_{i}.$

This integral is intractible, but you can compute the probability of interest numerically using the following Monte Carlo algorithm.

For $j$ from $1$ to $J$,

1. Sample a set of $z_i$ values from their posterior distribution.

2. Use the sampled values to compute $y_j=I(Q_{7}<Q_{6})I(Q_{7}<Q_{8})$ where $I(\cdot)$ is the indicator function.

Then $\Pr(Q_7<Q_6 \text{ and }Q_7<Q_8|f_{0},...,f_{8})\approx \frac{\sum_{j=0}^Jy_j}{J}$.

The accuracy of the Monte Carlo approximation goes as $\sqrt{J}$: $J=10^4$ will get you at least two decimal places of accuracy 19 times out of 20, $J=10^6$ will get you at least three decimal places of accuracy 19 times out of 20, etc.

And if your posterior probability of interest isn't close to 0 or 1, just sample more data, rinse, and repeat.

Prior hyperparameters, part two

The exponent of $z_i$ in the expression for the posterior density is

$\alpha_i + Nf_i - 1 = Ar^{i}\left(\frac{1-r}{1-r^{9}}\right) +Nf_i - 1 = A\mathrm{E}(z_i) +Nf_i - 1$

It can be seen that the hyperparameter $A$ plays the same role in the prior distribution as $N$ plays in the likelihood -- it's a kind of "prior sample size". To ensure that the prior has a negligible influence on the conclusion, just pick a value of $A$ such that $A\ll N$; for example, $A = 1$.

To set $r$, note that you can calculate the prior probability of the proposition $Q_7<Q_6 \text{ and } Q_7<Q_8$ using the same Monte Carlo algorithm described above but with the prior distribution in place of the posterior distribution in step 1 of the loop. Try to find a value of $r$ that gives a prior probability of 0.5 (or lower, if you feel that's more reasonable).

Answer 2

I presume the z_k are not probabilities but sample frequencies. This is because, otherwise, Q_i(z_0, ..., z_i) is not a random variable. In that case, computing the variance of the Q_i's is straightforward algebra. Define, first, the event indicators Z_i which is 1 if Z == i, 0 otherwise. It's a Bernoulli random variable with probability p_i. You can compute the first and second moments of any of these variables and they should give you all the necessary terms for computing the variance of the Q_i's.

Answer 3

Kevin, please be carefull since I'll have to change your notation a little bit: your $z_i$'s are not my $z_i$'s.

I think the following Bayesian solution is worth a try. Cook a random parameter $\Lambda>0$ and let $Z_1,\dots,Z_n$ be conditionally i.i.d., given $\Lambda=\lambda$, with $Z_i\mid\Lambda = \lambda \sim \textrm{Poisson}(\lambda)$. Use the notation $Z=(Z_1,\dots,Z_n)$. You already have a sample $z=(z_1,\dots,z_n)$ of the $Z_i$'s, with $n=2^{28}$. Define the random variables $$\Theta_i = P\{Z_i=k\mid \Lambda\} = \frac{e^{-\Lambda}\Lambda^k }{k!} \, , $$ for $i\geq 0$ (if this is not clear, take a look). Now, in this formulation your quadratic forms $Q_i=Q_i(\Theta_0,\dots,\Theta_i) = Q_i(\Lambda)$ are functions of $\Lambda$. So, the $Q_i$'s are random and you want to determine the posterior probability $$ P\{Q_7<Q_6 \,\,\,\textrm{and}\,\,\, Q_7<Q_8\mid Z=z\} \, . \qquad (*) $$ With a prior $\Lambda\sim\textrm{Gamma}(a,b)$, using Bayes Theorem we have $$ \Lambda\mid Z=z \sim \, \textrm{Gamma}\left( a + \sum_{i=1}^n z_i, b + n\right) \, . $$ You calculate $(*)$ generating i.i.d. $\lambda_i$'s from the former distribution (use R!) and computing $$ \frac{1}{N} \sum_{i=1}^N I_{(-\infty,Q_6(\lambda_i))\cap(Q_8(\lambda_i),\infty)}(Q_7(\lambda_i)) \, , $$ which converges, by the strong law of large numbers, to $(*)$ almost surely. To get a "yes" to your original question, this posterior probability must be "big enough". With such a huge sample ($n=2^{28}$), I think it is possible to play with the values of $a$ and $b$ to make your prior choice not much "informative".