Fast Dirichlet Distribution Parameter Estimation

From the Bayesian perspective, looking at \(\operatorname{Dir}(\alpha_{1}, \dots, \alpha_{n})\) the parameters \(\alpha_{j}\) encode the amount of evidence or number of observations of rolling an n-sided die, we witness side j approximately \(\alpha_{j}\) times.

It is tempting to use the Dirichlet distribution for estimating voting behavior, because a voter casts their ballot for a Democrat, a Republican, or a third-party candidate. The problem with this is, if we just set \(\alpha_{1}\) to the number of votes for the Democratic candidate (and likewise for the other parameters), there is too little variation. If we've seen a die turn up 1 and 6 roughly 99% of the time, we can be confident it will do so in the future. (We have such confidence about the Sun rising, after all.)

The problem: given observed outcomes \(x_{j}\) (dice rolls), we can form the approximate probabilities \(p_{j} = x_{j}/\sum_{\ell}x_{\ell}\). Using these probabilities, we can estimate the parameters of the Dirichlet distribution. I'm really interested in finding sensible values for \[\alpha_{0} = \alpha_{1} + \dots + \alpha_{n}.\tag{1}\] This quantity in Eq (1) is called the Precision of the Dirichlet distribution.

More precisely, we will observe an election in K counties, so we will end up with a vector \(\vec{x}_{j}\) whose components are the votes received by candidate j in each county. This gives us a vector \(\vec{p}_{j}\) whose components are the proportion of the county's vote went to candidate j

Solution: we can make estimates of the Dirichlet distribution based on the data, then use this to estimate the various parameters of interest.

Different Algorithms

There are a variety of ways to estimate parameters, I'm looking at closed-form algorithms. In particular, the method of moments derived algorithms.

In a nutshell, the method of moments looks at the moments of random variables \(\mathbb{E}[X^{n}]\) in terms of the parameters of the distribution, then tries to solve for one of the parameters in terms of the expected values. When given a population, we substitute in the sample mean (and sample variance) for the expected values in these formulas, giving us numerical estimates.

Minka's Method

Minka (2012) gives the following scheme. The method of moments applied to the Dirichlet distribution with parameters \(\alpha_{i}\), if we let \(\alpha_{0} = \alpha_{1}+\dots+\alpha_{n}\), may be approximated from data \(\vec{p}_{i}\) (in the trials the proportion corresponding to \(\alpha_{i}\)) using \[\widehat{\alpha_{0}} = \frac{\mathbb{E}[p_{1}] - \mathbb{E}[p_{1}^{2}]}{\mathbb{E}[p_{1}^{2}] - \mathbb{E}[p_{1}]^{2}}\] we find the approximation \[\boxed{\widehat{\alpha_{i}} = \mathbb{E}[p_{i}]\widehat{\alpha_{0}} = \mathbb{E}[p_{i}]\frac{\mathbb{E}[p_{1}] - \mathbb{E}[p_{1}^{2}]}{\mathbb{E}[p_{1}^{2}] - \mathbb{E}[p_{1}]^{2}}}\tag{2}\] However, there is no reason why we should use \(p_{1}\) in determining \(\widehat{\alpha_{0}}\) instead of any other \(p_{i}\).

Note algebraically speaking: \[\mathbb{E}[X^{2}] = \operatorname{var}[X] + \mathbb{E}[X]^{2}\]

A simple R implementation:

minka_estimate <- function(data) {
    mu <- mean(data[,1]);
    s2 <- var(data[,1]);

    (mu*(1 - mu))/s2 - 1;
}

Ronning's Method

We can use all the data to estimate \(\alpha_{0}\), as Ronning (1989) suggested. We use \[\operatorname{var}[p_{j}] = \frac{\mathbb{E}[p_{j}](1 - \mathbb{E}[p_{j}])}{1 + \alpha_{0}}\] after some algebraic rearrangement, then taking the log of both sides \[\log(\alpha_{0}) = \log\left(\frac{\mathbb{E}[p_{j}](1 - \mathbb{E}[p_{j}])}{\operatorname{var}[p_{j}]}-1\right).\] We take the average of the right-hand side over j (the candidates) to get \[\boxed{\log(\alpha_{0}) = \frac{1}{n}\sum^{n}_{j=1}\log\left(\frac{\mathbb{E}[p_{j}](1 - \mathbb{E}[p_{j}])}{\operatorname{var}[p_{j}]}-1\right).}\tag{3}\] The astute reader will recognize this is the geometric mean of estimates.

As a sanity test that this makes sense, we recover the correct value for the Beta Distribution.

A simple R implementation:

ronning_estimate <- function(data) {
    assert_that(all(0 <= data) && all(data <= 1));
    assert_that(all(0.97 <= rowSum(data)));
    
    col_estimate <- function(col) {
        mu <- mean(col);

        if (mu == 0) return(0);
        
        log(mu*(1-mu)/var(col) - 1);
    }
    
    exp(mean(apply(data, 2, col_estimate)));
}

Crazy algebraic results

With 3 results (a "trinomial", if you will), we have an estimate I found by hand. Let \(\mu_{A}\), \(\mu_{B}\) and \(\mu_{C}\) be the sample means and \(\sigma^{2}_{B}\) the sample variance. Assuming \[ 0\lt \mu_{A}\lt 1,\quad\mbox{and}\quad 0\lt\mu_{B}\lt 1-\mu_{A}\] (and if this fails to hold, permute the labels "A", "B", "C" until it holds, if at all possible), then we can estimate \[\widehat{\alpha_{A}} = \left(\frac{(1 - \mu_{B})\mu_{B}}{\sigma^{2}_{B}} - 1\right)\mu_{A} \tag{4.1}\] \[\widehat{\alpha_{B}} = \widehat{\alpha_{A}}\frac{\mu_{B}}{\mu_{A}} \] \[\widehat{\alpha_{C}} = \widehat{\alpha_{B}}\frac{1 - \mu_{B}}{\mu_{B}} - \widehat{\alpha_{A}}\]

Observe this uses Minka's formula in estimating \(\widehat{\alpha_{0}}\) (the parenthetic term on the right hand side of Eq (4.1)) for use in computing \(\widehat{\alpha_{A}}\). The reader can verify if we add these quantities together, we in fact get Minka's estimate using the second largest sample mean (and its associated sample variance).

Testing Results

When testing which estimator to use, we should have some error implicit in our mind, in the sense of "What are we trying to optimize? What counts as 'wrong'?" It's one thing to try to recover the parameters of a Dirichlet distribution, it's another to estimate the concentration parameter for a heuristic.

We could note that the method of moments for the Beta distribution coincide with estimating \(\alpha_{0}\) for using Minka's method on any particular outcome (i.e., using \(\widehat{\alpha_{j}}\) for any particular, fixed j). It's not uncommon for Bayesians to use a Gamma distribution as a prior for the concentration parameter of a Beta distribution, which would give us some measure of its variance.

The method of moments is biased (slightly). Despite this, we can still use the Cramer-Rao bound to estimate the variance of our estimates.

For what its worth, when used on voting results, the method of moments produce similar estimates when using the proportion of votes the Republican and Democrat received. It's within the margin of error of the maximum likelihood estimates given by MASS::fitdistr, but is consistently smaller than maximum likelihood estimates. Since I'm advocating humility in making these estimates, I prefer the method of moments using the fancy-pants algebraic method I derived (Minka's method on the runner-up candidate's proportion of votes).

If I were less lazy, I'd include more detailed proofs about estimates and goodness-of-fit tests, but...meh, exercise for reader :)

References

1. Jonathan Huang, "Maximum Likelihood Estimation of Dirichlet Distribution Parameters". Unpublished manuscript(?), undated, Eprint
2. Thomas P Minka, "Estimating a Dirichlet distribution". Microsoft Research Paper, 2012. Eprint.
3. G. Ronning, "Maximum-likelihood estimation of dirichlet distributions". Journal of Statistical Computation and Simulation 32 (1989) 215–221.

Sample Size and Central Limit Theorem

I am trying to reproduce results from Zachary R. Smith and Craig S. Wells's Central Limit Theorem and Sample Size.

Therein the authors claim the central limit theorem does not hold for samples of size 30 or so (they go so far as to claim 300). I have tried reproducing this claim, their work is unreproducible. Moreover, their work is visually wrong, you can literally plot the sums of uniformly distributed variables and observe the sum is normally distributed.

For those interested in the sordid details, I have posted it on github.

Why did they get this so wrong? Well, they failed to adequately handle multicomparison problems, and they ignored the issues involved with the Kolmogorov-Smirnov test. The latter was particularly troublesome, as it led them to completely erroneous findings.

Statistics as Decision Problem

"Decision theory" is a framework for picking an action based on evidence and some "loss function" (intuitively, negative utility). Almost all of statistics may be framed as a decision theoretic problem, and I'd like to review that in this post.

(Note that the diagrams in this post were really inspired by I-Hsiang Wang's lectures on Information Theory.)

I am going to, literally, give a "big picture" of statistics as decision theory. Then I'll try to drill down on various "frequentist" statistical tasks, to show the algorithmic structure to each one. Although I'm certain Bayesian data analysis can be made to fit this mould, I don't have as compelling a "natural" fit as frequentist statistics.

And just to be clear, we "the statistician" are "the decider" in this decision making problem. We are applying decision theory to the process of "doing statistics".

Review of Decision Theory

Statistical Experiment

We have some source of data, whether it's observation, experiment, whatever. As Richard McElreath's Statistical Rethinking calls it, we work in the "small world" of our model, where we describe the data as a random variable \(X\) which follows a hypothesized probability distribution \(P_{\theta}\) where \(\theta\) is a vector of parameters describing the "state of the world" (it really just parametrizes our model). The set of all possible parameters is denoted \(\Theta\) with a capital theta. This \(\Theta\) is the boundary to our "small world". This data collection process is highlighted in the following figure:

Serious statisticians need to actually think about sampling methods and experimental methods. We are silly statisticians not serious statisticians, and use data already assembled for us. Although we will not perform any polling or statistical experiments ourselves, it is useful to know the nuances and subtleties surrounding the methodology used to produce data. Hence we may dedicate a bit of space to discuss the aspects of data gathering and experimental methodologies our sources have employed.

Decision Making

Given some data we have just collected, we now arrive at the romantic part of data collection: decision making. Well, that's what decision theorists call it. Statisticians call it...well, it depends on the task. It is highlighted in the following diagram:

There are really multiple tasks at hand here, so lets consider the key moments in decision making.

Inference task. Or, "What do we want to do with the data?" The answer gives us the task of estimating a specific function \(T\) of the parameters \(T(\theta)\) from the observed data X. The choice of this function depends on the task we are trying to accomplish.

A few examples:

• With hypothesis testing, \(T(\theta)=\theta\) we're trying to estimate the parameters themselves (which label the hypotheses we're testing).
• For regressions (i.e., given pairs of data from the experimental process \((X,Y)\) find the function f such that \(Y = f(X) + \varepsilon\)) the function of the parameters is the relationship itself, i.e., \(T(\theta)=f\).
• For classification problems, \(T(\theta)\) gives us the "correct" labeling function for the data.

In some sense, \(T(\theta)\) is the "correct answer" we're searching for, we just have to approximate it with the next step of the game...

The Decision Rule. In the language of decision theory, an estimator is an example of a Decision Rule which we denote by \(\tau\) ("tau"). This approximates \(T(\theta)\) given the data we have and the conceptual models we're using.

For regressions, this is the estimated function \(\tau(X,Y)=\widehat{f}_{X,Y}\) which fits the observations. For hypothesis testing, \(\tau(X)=\hat{\theta}\) is which hypothesis "appears to work".

These two tasks, inference and computing the decision rule, constitutes the bulk of statistical work. But there's one more crucial step to be done.

Performance Evaluation

We need to see how good our estimates are! In the complete diagram, this is the highlighted part of the following figure:

The loss function \(l(T(\theta),\tau(X))\) measures how bad, given the data X, the decision rule \(\tau\) is. Note this is a random variable, since it's a function of the random variable X. Also note, there are various different candidates for the loss function (it's our job as the statistician to figure out which one to use).

The risk is just the expected value of the loss function. This tells us on average how bad the decision rule \(\tau\) turns out, given the true state of the world is \(\theta\). We denote this risk by the function \(L_{\theta}(\tau)\).

For some tasks, we don't really have much of a choice on the loss function. Regressions do best with the mean-squared error. We could choose a different loss function (e.g., a variant on the mean squared error, we could use the \(L^{p}\) norm instead of the \(L^{2}\) norm).

Remark. It might seem strange that, given \(T\) is never really knowable, yet it appears in the risk function. We typically use it in a tricky way. For example in regression we're really using \(T(\theta)=f\) and then the trick is use \(f(X) = Y\) and we can use the observed results Y instead of worrying about the unobservable, incalculable, unknowable \(T\).

Examples

We will collect a bunch of examples, but this is incomplete. The goal is to show enough examples to encourage the reader to devise their own.

Hypothesis Testing

Classical hypothesis testing may be framed as a decision problem: do we take action A or action B? For our case, do we accept or reject the null hypothesis.

More precisely, we have two hypotheses regarding the observation X, indexed by \(\theta=0\) or \(\theta=1\). The null hypothesis is that \(X\sim P_{0}\), while the alternative hypothesis states \(X\sim P_{1}\).

We have some decision rule, which in our diagrams we have denoted \(\tau(X)\), which "picks" a \(\theta\) which minimizes the risk based on the observations X. But what is the loss function?

Well, we have the probability for a false alarm when \(\tau(x)=1\) but should be zero \[\alpha_{\tau} = \sum_{x}\tau(x)P_{0}(x)\tag{1}\] and the probability for missing a detection when \(\tau(x)=0\) but should be one \[\beta_{\tau} = \sum_{x}(1-\tau(x))P_{1}(x)\tag{2}.\] We note the loss function is indeed an expected value of \(\tau\), and it is parametrized by the choice of \(\theta\).

But how do we choose \(\tau\)?

We may construct one possible "hypothesis chooser" (randomized decision rule) as, for some constant probability \(0\leq q\leq 1\) and threshold \(c\gt0\), \[\tau_{c,q}(x) = \begin{cases} 1 & \mbox{if } P_{1}(x) \gt cP_{0}(x)\\ q & \mbox{if } P_{1}(x) = cP_{0}(x)\\ 0 & \mbox{if } P_{1}(x) \lt cP_{0}(x) \end{cases}\tag{3}\] In other words, \(\theta=1\) is chosen with probability \(\tau_{c,q}(x)\), and \(\theta=0\) is chosen with probability \(1-\tau_{c,q}(x)\). Starting from a given value of \(\alpha_{0}\), we then determine the parameters c and q by the equation \[\alpha_{0}=\sum_{x}\tau_{c,q}(x)P_{0}(x).\tag{4}\] The Neyman-Pearson lemma proves this is the most powerful test for significance (minimizes \(\beta_{\tau_{c,q}}\) subject to \(\alpha_{\tau_{c,q}}=\alpha_{0}\) constraint).

We emphasize, though, this is a "toy problem" which fleshes out the details of this framwork.

Exercise. Prove that the probability of type-I errors (probability of false alarm) is \(\alpha_{\tau} = \mathbb{E}_{X\sim P_{0}}[\tau(X)]\) and the probability of type-II errors (probability of failing to detect) is \(\beta_{\tau} = \mathbb{E}_{X\sim P_{1}}[1 - \tau(X)]\).

Regression

The goal of a regression is, when we have some training data \((\mathbf{X}^{(j)}, Y^{(j)})\) where parenthetic superscripts run through the number of observations \(j=1,\dots,N\), to find some function f such that \(\mathbb{E}[Y|\mathbf{X}]\approx f(\mathbf{X},\beta) \approx Y\). Usually we start with some preconception like f is a linear function, or a logistic function, or something similar, rather than permitting f to be any arbitrary function. We then proceed to estimate \(\widehat{f}\) and the coefficients \(\widehat{\beta}\).

Some terminology: the \(\mathbf{X}\) are the Covariates (or "features", "independent variables", or most intuitively "input variables") and \(Y\) are the Regressands (or "dependent variables", "response variable", "criterion", "predicted variable", "measured variable", "explained variable", "experimental variable", "responding variable", "outcome variable", "output variable" or "label"). Unfortunately there is a preponderance of nouns for the same concepts.

Definition. Consider \(X\sim P_{\theta}\) which randomly generates observed data \(x\), where \(\theta\in\Theta\) is an unknown parameter. An Estimator of \(\theta\) based on observed \(x\) is a mapping \(\phi\colon\mathcal{X}\to\Theta\), \(x\mapsto\hat{\theta}\). An Estimator of a function \(z(\theta)\) is a mapping \(\zeta\colon\mathcal{X}\to z(\Theta)\), \(x\mapsto\widehat{z}\).

The decision rule then estimates the true function, \(\tau_{\mathbf{X},Y}=\widehat{f}\). That is to say, it produces an estimator. There are various algorithms to construct the estimator, which depends on the regression analysis being done.

The loss function is usually the squared error for a single observation \[l(T,\tau) \mapsto \mathbb{E}_{(\mathbf{X},Y)\sim P_{\theta}}[(Y - \widehat{f}(\mathbf{X},\widehat{\beta}))^{2}].\tag{5}\] Depending on the problem at hand, other loss functions may be considered. (If the Y variables were probabilities or indicator variables, we could use the [expected] entropy as the loss function.)

The risk is then the average loss function over the training data. But do not mistake this for the only diagnostic for regression analysis.

We have multiple measures of how good our estimator is, which we should briefly review.

Definition. For an estimator \(\phi(x)\) of \(\theta\),

• its Bias is  \(\mathrm{Bias}_{\theta}(\phi) := \mathbb{E}_{X\sim P_{\theta}}[\phi(X)] - \theta\)
• its Mean Square Error is  \(\mathrm{MSE}_{\theta}(\phi) := [|\phi(X) - \theta|^{2}]\)

Fact (The MSE = (Bias)² + Variance). Let \(\phi(x)\) be an estimator of \(\theta\), then \[\mathrm{MSE}_{\theta}(\phi) = \left(\mathrm{Bias}_{\theta}(\phi)\right)^{2} + \mathrm{Var}_{P_{\theta}}[\phi(X)]. \tag{6}\] In practice, this means as an estimate is more "spread out", it becomes more "centered near the correct value". (End of Fact)

Conclusion

Most of statistical inference falls into this schema presented. Broadly speaking, statistical inference consists of hypothesis testing (already discussed), point estimation (and interval estimation), and confidence sets.1See, e.g., section 6 of K.M. Zuev's lecture notes on statistical inference. We have discussed only the frequentist approach, however, and for only a couple of these tasks.

The Bayesian approach, in contrast to all these techniques, end up using a loss function which sums over values of \(\theta\) (i.e., integrates over \(\Theta\) instead of over the space of experimental results \(\mathcal{X}\)). The Bayesian priors describe a probability distribution of likely values of \(\theta\), which would be used in the overall process.

Yet the Bayesian school offers more tools than just this, and I don't think they can neatly fit inside a diagram like the one doodled above for frequentist statistics.

But we have provided an intuition to the overall procedures and tools the frequentist school affords us. Although we abstracted away the data gathering procedure (as well as the other steps in the process), we could flesh out more on each step involved.

In short, statistics consists of decision theoretic problems (perhaps "decision theory about decision theory", or meta-decision theory, may be a good intuition), but it remains more of an art than an algorithmic task.

References

• P.J.Bickel and E.L.Lehmann, "Frequentist Inference". In Neil J. Smelser and Paul B. Baltes (eds) International Encyclopedia of the Social & Behavioral Sciences, Elsevier, 2001, Pages 5789–5796.
• James O. Berger, Statistical Decision Theory and Bayesian Analysis. Springer Verlag, 1993. See especially section 2.4.3. (This is the only book on statistical decision theory that I know of worth its salt.)
• George Casella and Roger L. Berger, Statistical Inference. Second ed., Cengage Learning, 2001. (Section 8.3.5, 9.3.4 for hypothesis testing and point estimators in the decision theoretic framework.)
• C. Robert, The Bayesian choice: from decision-theoretic foundations to computational implementation. Second ed., Springer Verlag, 2007. See especially chapter 2.

Exit Polls Margin of Error Estimates

Exit polls do not have margins of error, but we can estimate the margin of error using confidence intervals for a two-party system.

Applied to the 2016 election, when we are told N respondents answered with a proportion p supporting a candidate, we can construct the Wilson confidence interval with, say, 95% confidence (i.e., α = 0.05 and \(z_{1-\alpha/2}\approx 1.96\)). Then we get an estimate \(\hat{p}\pm\Delta p\), and we may treat \(\Delta p\) as the margin of error.

If we are trying to estimate uncertainty propagated from the exit polls used in, say, computing coefficients for a logistic regression, then we could use z = 2 exactly, and then set the Wilson confidence interval computed to \(\hat{p}_{A}\pm2\sigma_{A}\) for supporters of candidate A.

The only caveat is, exit polls are not adequately random samples. Exit polls are cluster samples, since only a fraction of precincts are polled (although they are picked by random and are intended to reflect the state as a whole). There are techniques for computing the margin of error for such sampling techniques, I don't believe it to be tractable given the limited data from exit polls.

We can compute the margin of error for one-stage cluster sampling (which would be the upper bound in the margin of error, i.e., the stratified cluster sampling would have a smaller margin of error). How does it compare to binomial confidence intervals? Lets review binomial confidence intervals, then cluster sampling error, and see when/if the binomial confidence interval is a superior choice.

Brief Review of Confidence Intervals for Binomial Distribution

Remember the de Moivre-Laplace theorem, which states if X is a binomially distributed random variable with probability p of success in n trials, then as \(n\to\infty\) we find \[\frac{X-np}{\sqrt{np(1-p)}}=Z\to\mathcal{N}(0,1) \tag{1} \] the left hand side becomes approximately a normal distribution with mean 0 and standard deviation 1. Dividing top and bottom by n, then rearranging terms, we find \[\frac{X}{n} = p + Z\sqrt{p(1-p)/n}\tag{2}\] which gives us an estimate for p.

If we denote \(\hat{p}=y/n\) the empirically observed frequency of successes (y) to the number of trials (n), we pick some confidence level z, and the naive interval estimate for the probability of success is given by the normal approximation \[\boxed{\widehat{p}\pm z\sqrt{\frac{\widehat{p}(1-\widehat{p})}{n}}}\tag{3}\] which, for large enough n and for \(\hat{p}\) "not too extreme", gives some estimate for where the "true value" of p lies.

Puzzle: What values of n and p lead to good approximations by Eq (3)?

This puzzle is typically "solved" in textbooks by insisting \(n\cdot\mathrm{min}(\hat{p},1-\hat{p})>5\) (or 10), but this doesn't always lead to good estimates. Brown, Cai, and DasGupta investigated this question "empirically".

In fact, a better estimate of the interval starts by considering p at the boundaries of Eq (3) with a slight modification of the squareroot: \[p = \widehat{p}\pm z\sqrt{\frac{p(1-p)}{n}}\tag{4}\] then we get the quadratic equation \[(p - \widehat{p})^{2} = z^{2}\frac{p(1-p)}{n}.\tag{5}\] Solving the quadratic for p gives us the interval \[\boxed{\frac{\hat p+\frac{z^2}{2n}}{1+\frac{z^2}{n}} \pm \frac{z}{1+\frac{z^2}{n}}\sqrt{\frac{\hat p(1-\hat p)}{n}+\frac{z^2}{4n^2}}.}\tag{6}\] This is the Wilson Confidence Interval. Heuristically, for \(z\approx 2\) (the 95% confidence interval) estimates the true probability of success (p) to be centered nearly at a shifted estimate \((y+2)/(n+4)\). Observe for "large n", Eq (6) becomes Eq (3).

The Wilson confidence interval gives better estimates than the normal approximation, even for a small number of trials n and/or extreme probabilities. For larger n (i.e., \(n\gt 40\)), the Agresti-Coull interval should be used: compute the center of the Wilson confidence interval, then use this value as \(\hat{p}\) in the normal approximation Eq (3).

In short: If \(n\lt 40\), use the Wilson confidence interval. Otherwise, use either the Agresti-Coull interval, the Wilson interval, or the usual interval Eq (3).

Remark. For small n < 40, we could use a Bayesian approximation with the uninformative Jeffreys prior. This estimates the interval, for a confidence level α, to be the quantiles of the Beta distribution \(\mathrm{Beta}(y + 1/2, n-y+1/2)\) at the probabilities \(\alpha/2\) and \(1-\alpha/2\). This has to be computed numerically. I wonder if there are decent approximations to the quantile function for small α?

Rare Events

What if the probability of success is really low? That is to say, we're dealing with "rare events"? We have a special case for this, going back to the binomial distribution to describe the most extreme case where there are no successes observed y = 0 is described by \[\Pr(X=0) = (1 - p)^{n} = \alpha\tag{7}\] for a given confidence level α (usually 0.05). Then taking the (natural) logarithm of both sides yields \[n\ln(1-p)=\ln(\alpha).\tag{8}\] By assumption, the chance for success is small (\(p\ll 1\)), we can approximate the logarithm by the first term of the Taylor expansion (since the first term is an upper bound of the logarithm in this domain) \[-np=\ln(\alpha)\tag{9}\] hence the confidence interval is \[\boxed{0\leq p\leq\frac{-\ln(\alpha)}{n}.}\tag{10}\] For α = 0.05, the upper bound is 3/n (hence the so-called "Rule of three").

(Dually, for events which almost always happen, we could simply take 1 minus this interval. For α = 0.05, this is nearly [1-3/n,1].)

Cluster Sampling Margin of Error

The cluster sampling margin of error is the product of the standard error with the critical value z. We are given the sample proportions p. The cluster, in the case of exit polling, is a precinct; in 2004, there were 174,252 precincts in the United States (arXiv:1410.8868) with an average of 800 voters in a precinct. There is something on the order of 300 precincts (in 28 states, roughly 11 per state) sampled in the 2016 exit poll.

If there are N voters in the country, N_j voters in precinct j, m precincts in the exit poll, and M precincts in the country, and y_j be the number of voters in precinct j that voted for a fixed party for president, then we may consider the estimated number of votes may be given using the unbiased estimator \[\hat{\tau} = M\cdot\bar{y} = \frac{M\cdot\sum^{m}_{j=1}y_{j}}{m}.\tag{11}\] Its variance is given by \[\mathrm{Var}(\hat{\tau}) = M(M-m)\frac{s_{u}^{2}}{m}\tag{12a}\] where \[s_{u}^{2} = \frac{1}{m-1}\sum^{m}_{j=1}(y_{j}-\bar{y})^{2} \tag{12b}\] is the sample variance.

We can estimate the proportion of voters supporting our given party. We take \[\hat{\tau}_{r} = N\cdot r = N\cdot\frac{\sum^{m}_{j=1}y_{j}}{\sum^{m}_{j=1}N_{j}}\tag{13a}\] and \[\hat{\mu}_{r} = \hat{\tau}_{r}/N = r.\tag{13b}\] We find the variance \[\mathrm{Var}(\hat{\tau}_{r}) = \frac{M(M-m)}{m(m-1)}\sum^{m}_{j=1}(y_{j} - rN_{j})^{2}\tag{14a}\] which is biased, but the bias is small when the sample size is large; the variance for the ratio estimator \[\mathrm{Var}(\hat{\mu}_{r}) = \frac{M(M-m)}{m(m-1)}\frac{1}{N^{2}}\sum^{m}_{j=1}(y_{j} - rN_{j})^{2}.\tag{14b}\] It is not hard to find \[\mathrm{Var}(\hat{\mu}_{r}) = \frac{(1-m/M)}{m(m-1)}\sum^{m}_{j=1}\left(\frac{y_{j}}{N_{j}} - r\right)^{2}\frac{N_{j}^{2}M^{2}}{N^{2}}.\tag{14c}\] which is the variance we are looking for.

The margin of error for the exit polls would be approximately, for a given critical value z, \[ME = z\sqrt{\mathrm{Var}(\hat{\mu}_{r})}.\tag{15}\] Unfortunately, we are not given the data sufficient to compute the variance described in Eq (14c).

Another difficulty, Eq (14c) describes the sample variance. Exit polls are far less than ideal (in the US, at least), and this increases the actual variance. There has been some debate surrounding how much worse the error for exit polls is, when compared to naive binomial confidence interval estimates, but the Mystery Pollster inform us, Panagakis had checked with Warren Mitofsky, director of the NEP exit poll, and learned that the updated design effect used in 2004 assumed a 50% to 80% increase in error over simple random sampling (with the range depending on the number of precincts sampled in a given state). (Emphasis his) If the reader takes one thing away from this post, it should be exit polls are noisy and computing its margin of error is complicated.

Conclusion: multiplying the width of the confidence interval by a factor of 1.8, or even 2, would give us a reasonable margin of error for the exit polls.

References

Stat506 from Pennsylvania State University is where I learned about cluster sampling.

2016 Exit Polls

Exit polls, by their nature, are extremely noisy and do not provide margins of error. Until 2016, news organizations banded together to provide a single, coherent exit poll for presidential elections. This coalition started collapsing after the 2016 presidential election. Both CNN and Fox News have slightly different exit poll data for the 2016 election. How do we make sense of these polls? When can we say they tell us "the same story"?

Toy Problem: Coin Tossing

Lets consider a simpler toy problem: I flip a coin N₁ times and obtain y₁ heads ("successes") and you flip a coin N₂ times and obtain y₂ heads ("successes"). How do we know our coins are equally biased? Assuming each of us has done a "large number of trials" (a couple hundred each).

The null hypothesis: these coins follow the same distribution with probability p of success.

The alternative hypothesis: these coins do not follow the same distribution.

Let \(p_{1} = y_{1}/N_{1}\) (and similarly \(p_{2} = y_{2}/N_{2}\)), a better approximation to the probability of heads ("success") is \[p=\frac{y_{1}+y_{2}}{N_{1}+N_{2}}.\tag{1}\] We can standardize the data and approximate it as a normal distribution with standard deviation \(\sigma_{1}^{2} = p(1-p)/N_{1}\) and similarly \(\sigma_{2}^{2} = p(1-p)/N_{2}\). Then the "true variance" is \[\sigma^{2} = \sigma_{1}^{2} + \sigma_{2}^{2} = p(1-p)\left(\frac{1}{N_{1}}+\frac{1}{N_{2}}\right)\tag{2}\] which we should use for performing the Z-transform: \[Z_{1} = \frac{p_{1} - p}{\sigma}\tag{3}\] and a similar definition of Z₂. The test statistic we want to measure is \[z = Z_{1}-Z_{2} = \frac{p_{1}-p_{2}}{\sigma}\tag{4}\] which should follow a normal distribution with mean 0 and standard deviation 1.

We can then pick some confidence level, usually 1 - α = 95% confidence level, which leads to the critical value \[z_{1 - \alpha/2} = \Phi^{-1}(1 - \alpha/2)\tag{5}\] using the quantile for the normal distribution \(\Phi^{-1}\) for approximately 1.96 for the 95% confidence level.

If the quantity computed in Eq (4) is "more extreme" than the quantity computed in Eq (5), i.e., if \(z_{1-\alpha/2} < |z|\), then we reject the null hypothesis and conclude these coins follow different distributions. Otherwise, we fail to reject the null hypothesis. Note: we never prove they follow the same distribution, we just fail to prove they follow different distributions.

Exit Polls

We can apply this same setup to exit polls. Unfortunately for Fox News, contrary to appearances, they do not provide state-level exit poll data for the first 12 questions. Perhaps this is a technical error due to some software bug on their server-side, but I can't do anything to remedy the problem.

Further, the results are identical percentages (on the first dozen questions) for FOX and CNN, despite having apparently different sample sizes. It is not hard to prove (if these results are accurately reported) the results of computing Eq (4) for each category is 0 identically, and thus we fail to reject the null hypothesis for each question.

It's rather anti-climactic despite this long buildup, but it's good to know these exit polls are coherent, in some appropriate sense.

The exit polls are in CSV form on github.

What is...a Random Variable?

Tentative Definition. A random variable assigns to a given "random phenomenon" or "random process" some (real) number, or a vector ("list") of numbers.

Examples. The following are all random variables.

1. Let X be the number of heads in 10 coin flips.
2. Let R be the number of times a given baseball pitcher strikes out an opponent in the course of a given game.
3. Let S be the number of "successes" (heads) until the first "failure" (tail) in a repeated trial (flipping a coin over and over again).
4. Let T be the waiting time (in minutes) until the next bus arrives.

Slightly more formal Definition. If we represent the possible outcomes for a given "random process" (i.e., the set underlying its sigma algebra) by Ω, then a random variable is a function \(X\colon\Omega\to\mathbb{R}\) (or possible \(X\colon\Omega\to\mathbb{R}^{k}\) some [fixed positive integer] k) such that for any \(x\in\mathbb{R}\) the preimage of smaller values is an event \(\{\omega\in\Omega:X(\omega)\leq x\}\in\Sigma\) ("is measurable", i.e., we can assign a probability to that preimage). We will write \(X\leq x = \{\omega\in\Omega:X(\omega)\leq x\}\) in an abuse of notation.

Remark. We should remember a sigma algebra is not just the space of outcomes \(\Omega\), but also the specific set of well-defined events \(\Sigma\subseteq\mathcal{P}(\Omega)\). There are various specifications we have on the events, for instance "something happens" \(\Omega\in\Sigma\); for any event \(E\in\Sigma\) its complement is also a well-defined event \(\Omega\setminus E\in\Sigma\) [or "an event does not happen"]; for any countable family of events \(\{E_{j}\}_{j\in J}\subseteq\Sigma\) their union is also a well-defined event \(\bigcup_{j\in J}E_{j} \in\Sigma\) ["one of these events might happen"]; and so forth. Implicitly we consider a probability measure on the set of well-defined events generically denoted \(\Pr(-)\). Altogether, this is the data necessary to describe some "random phenomenon".

Example. Let E be any event. Then indicator function \(I_{E}\) which is zero for any \(x\notin E\) and \(I_{E}(e)=1\) for all \(e\in E\), this is a random variable. For multiple events \(E_{1},\dots, E_{n}\), we find \(I_{E_{1}\cup\dots\cup E_{n}}(x)=\max(I_{E_{1}}(x),\dots,I_{E_{n}}(x))\) and \(I_{E_{1}\cap\dots\cap E_{n}}(x) = \prod_{j} I_{E_{j}}(x)\). For discrete random phenomena, these indicator functions are the basic building blocks for constructing other random variables.

What happens to all that sigma-algebra baggage? Given a random variable, we can ask "What is the probability its value will be in a given range?" For example, "What is the probability the starting pitcher for the Dodgers will strike out at least 30 batters?"

This would be computed by first assembling all possible outcomes which satisfy this \(\mathcal{E}=\{E\in\Omega : R(E)\geq30\}\subseteq\Omega\) and using the probability measure on the sample space for the random process \(\Pr(\mathcal{E})\) or more imaginatively \(\Pr(R\geq30)\).

As a caveat, we hasten to add, only inequalities are necessarily well-defined, e.g., \(\Pr(X\leq x)\). Equality "on the nose" may not be well-defined \(\Pr(X=x)\), but we abuse notation to write the Probability mass function in this manner. (This gets tricky with subtle nuances when dealing with continuous random variables instead of discrete ones.)

This induces a nice mathematical structure on the image of the random variable \(R(\Omega)\), namely we can "transport" the probability distribution from the sigma algebra \(\Sigma\) on \(\Omega\) to \(R(\Omega)\). This is the Distribution Function for the random variable, \(F_{R}(x) = \Pr(R\leq x)\).

Equivalence relation of Random Variables. If we have two random variables, say, \(X\) and \(Y\), we can say they are Equivalent if for any \(x\in\mathbb{R}\), we have \(\Pr(X\leq x) = \Pr(Y\leq x)\). This is usually denoted \(X\sim Y\).

Algebra of Random Variables. Given some random variables X, Y on the sample sigma algebra, we can define new random variables \(X+Y\), \(X - Y\), \(XY\), \(X/Y\) provided Y is never zero, and exponentiation \(X^{Y} = \exp(Y\log(X))\). The intuition to have is that the operations are done as operations on real-valued functions.

So, specifically, if \(x\in\Omega\), then \((X+Y)(x)=X(x)+Y(x)\), \((X - Y)(x) = X(x) - Y(x)\), \((XY)(x) = X(x)Y(x)\), \((X/Y)(x) = X(x)/Y(x)\) provided \(Y(x)\neq0\), and exponentiation \((X^{Y})(x) = \exp(Y(x)\log(X(x)))\).

Probability Distributions. We have a few "standard" distributions which are the "template" for various random processes. Flipping a coin follows a Bernoulli Distribution if we flip the coin only once, and a Binomial Distribution if we flip it N times, for example.

The notation for these families may vary reference to reference. A Bernoulli distribution with probability p of success is usually denoted \(\mathrm{Bernoulli}(p)\) or \(\mathrm{Ber}(p)\).

To indicate a random variable is distributed like one of these standard distributions, we abuse notation and write \(X\sim \mathrm{Bernoulli}(p)\).

We can build more distributions out of a handful of basic ones, for example \(Y = X_{1} + \dots + X_{N}\) where the \(X_{j}\sim\mathrm{Bernoulli}(p)\) will describe flipping a coin N times and counting the number of "successes" ("heads"). This gives us the Binomial distribution when we consider \(\Pr(Y\leq k)\) (there are at most k heads in N coin flips).

We can specify a probability distribution by its parameters (like p in the Bernoulli distribution), the probability mass function and/or the probability density function. Often it's useful to give other summary statistics alongside this data.

Expected Value. We also have for any random variable X its Expected Value given by \(\mathbb{E}[X] = \sum_{x\in X(\Omega)}x\Pr(X=x)\) (or replacing the sum with an integral for continuous random variables). The intuition we should have for the expected value of a random variable is this captures the "average value" of the random variable.

If we have some function \(f\colon\mathbb{R}\to\mathbb{R}\), then we have \(\mathbb{E}[f(X)] = \sum_{x\in X(\Omega)}f(x)\Pr(X=x)\) (and again, an integral instead of a sum for continuous random variables, with the restriction that f is an integrable function).

Note that \(\mathbb{E}[X^{2}]\neq(\mathbb{E}[X])^{2}\) and more generally \(\mathbb{E}[XY]\neq \mathbb{E}[X]\mathbb{E}[Y]\). But we do have \(\mathbb{E}[X + Y] = \mathbb{E}[X] + \mathbb{E}[Y]\) and, for any real number \(a\in\mathbb{R}\) \(\mathbb{E}[aX] = a\mathbb{E}[X]\).

Exercise. Let \(E\in\Sigma\) be an event in a sigma algebra, and \(I_{E}\) be the indicator function on E. What is \(\mathbb{E}[I_{E}]\)?

Theorem. Let X and Y be random variables, a and b be real numbers. Then:

1. \(\mathbb{E}[aX + b] = a\mathbb{E}[X] + b\)
2. \(\mathbb{E}[X+Y] = \mathbb{E}[X] + \mathbb{E}[Y]\)
3. \(\displaystyle\mathbb{E}[XY] = \sum_{\omega\in\Omega}X(\omega)Y(\omega)\Pr(\omega)\)
4. \(\displaystyle\mathbb{E}[X/Y] = \sum_{\omega\in\Omega}\frac{X(\omega)}{Y(\omega)}\Pr(\omega)\) if  \(Y(\omega)\neq0\) for any \(\omega\in\Omega\)
5. \(\displaystyle\mathbb{E}[X^{Y}] = \sum_{\omega\in\Omega}\exp(Y(\omega)\log(X(\omega)))\Pr(\omega)\)

Variance. If expected value tells us what neighborhood a random variable is likely to live in, the variance tells us how spread out this neighborhood is. We define it as \[\mathrm{Var}[X] = \mathbb{E}[X^2] - (\mathbb{E}[X])^2 = \sum_{\omega}(X(\omega) - \mathbb{E}[X])^{2}\Pr(\omega).\] The variance and expected value for a random variable contain a lot of useful information, which we use when trying to infer parameters from data.

More generally, we have for any two random variables X and Y a measure of their failure to be correlated by the covariance \[\mathrm{Cov}(X,Y) = \sum_{\omega}(X(\omega)-\mathbb{E}[X])(Y(\omega)-\mathbb{E}[Y])\Pr(\omega)\] which is such that \(\mathrm{Cov}(X,X)=\mathrm{Var}[X]\). Correlatedness is not the same as independence: independent random variables are not correlated, but uncorrelated random variables may or may not be independent (a "all salmon are fish, not all fish are salmon" type statement). So uncorrelated is a "weaker" property than independence.

Applications?

Iterate! Note, we can transform the parameters of these distributions (like p the probability of success in a Bernoulli trial) into random variables themselves. This is precisely what Bayesian data analysis does: the parameters are random variables following prior probability distributions, which we update as new data becomes available using Bayes's theorem.

Regressions! A linear regression basically says that the observations are really values of a random variable, i.e., \(Y\sim\mathcal{N}(aX + b, \varepsilon)\) where \(\mathcal{N}\) is the normal distribution. There are other useful regressions, but this is the basic idea.

We should admit that this is one formulation of regressions in terms of random variables. The other uses conditional random variables, \(Y|X\sim f(\beta\cdot X,\theta)\) when the regressions (X) are stochastic (i.e., "not controlled by the experimenter/statistician"). Formally these are different models. But when actually doing the regressions, they are treated "the same".

Tests! We often do an experiment, producing some data points \(x_{1},\dots,x_{n}\) which we interpret as values of a random variable X which follows a prescribed distribution. We test the assumption (that X follows the given distribution with specific parameters) by comparing the sample mean \[ \mu = \frac{1}{n}\sum_{j=1}^{n}x_{j} \] to the expected value \(\mathbb{E}[X]\). The central limit theorem suggests that \(\sqrt{n}(\mu - \mathbb{E}[X])\) looks like a normal distribution centered at 0 with variance approximately equal to the variance of the data points (loosely speaking).

Reading List

• Discrete Random Variables and their Expected Values, MIT OCW 18-05 Introduction to Probability and Statistics
• Expected Value and Variance, lecture notes.
• Statistical Machine Learning, by Han Liu and Larry Wasserman, specifically chapter 12 on Bayesian inference using random variables.