# Confidence intervals for a polynomial

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.

• Your uneasiness with the language of confidence intervals, but otherwise rigorous presentation, suggests you may be in a good position to give us a little more information to help you. Ordinarily you can specify one of two properties of a CI: either its coverage or its length. The other will be determined by the data. It's unusual to stipulate the length, though: are you sure this is what you want? Also, it is unclear whether you need CIs for the Q's separately or simultaneously. Perhaps you could indicate what you intend to use these CI's for? – whuber Feb 27 '12 at 23:32

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

• Naw, I'm just skipping some technicalities. A more technically complete analysis would start with a Dirichlet process and then show that the result of marginalizing out the countably infinite set of irrelevant $z_i$ parameters is the Dirichlet distribution I give above. – Cyan Mar 30 '12 at 2:57

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.

• No, the $z_i$ are constants of nature, and consequently $Q_7$ is, too. It's a real number, but the question is which real number. My sampling lets me estimate the $\hat{z}_i$, and if I plug those into the formula linking $Q_7$ to the $z_i$, I get an estimate for $Q_7$ (around $0.07$). But I don't understand how good that estimate is. If I take 99% CIs for each $z_i$ and plug those in, I get an interval, but what is the confidence level? I don't think the CIs for the $z_i$ are independent, either. – Kevin O'Bryant Feb 27 '12 at 23:07
• Lack of independence isn't a worry here, Kevin (you can justify this by analyzing the true--multinomial--distribution of the $\hat{z}_i$ if you wish). – whuber Feb 27 '12 at 23:29
• Since you can already compute an estimate for $Q_7$ from your data, how about using some resampling technique like bootstrap to find a CI for $Q_7$? en.wikipedia.org/wiki/Bootstrapping_(statistics) – Zen Feb 28 '12 at 2:28
• @Zen: I'd heard of bootstrapping, but never thought it would be something that happened to me. Seriously, I'm trying to make it work, but it's time consuming to resample 1000 times (each resampling having $2^{28}$ points). Even 1000 resamples of size $2^{16}$ takes me 2 hours. – Kevin O'Bryant Mar 13 '12 at 15:36

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".

• That Poisson distribution seems like a rather restrictive assumption, no? – Cyan Mar 29 '12 at 23:27
• We need a distribution, whose support is $\mathbb{Z}_+$, that makes the computations possible. Let's see what Kevin can find with this model. – Zen Mar 30 '12 at 2:14