1. Make a table for the joint distribution of three rolls of a 3-sided die. Add columns for the responses.
2. Partially marginalize to get the joint distribution of the responses. Add a column for their sum.
3. Marginalize to get the distribution of the sum.

• To find the distribution of a sum of two responses, \(Y_1+Y_2\), we do the same thing.

• We start with the joint distribution of two dice rolls.

  • Our rolls are equally likely to be any pair of numbers in \(1\ldots m\).
  • And there are \(m^2\) pairs, so the probability is \(1/m^2\) for each pair.

• Then we add columns that are functions of the rolls: \(Y_1\), \(Y_2\), and \(Y_1+Y_2\).

  • This is ‘for free’ as before.
  • If we know who we’re calling, we know the responses we’ll hear and their sum.

• Then we marginalize to find the distribution of the sum. We’ll do this in two steps.

  1. We partially marginalize over \(J_1,J_2\), to find the joint distribution of the response sequence \(Y_1, Y_2\).
  2. We fully marginalize the result to find the distribution of \(Y_1+Y_2\).

• To sum over pairs of responses, we sum over one response then the other.
  • This exactly like what we did for one response. Except twice.
  • And what we get is pretty similar too.
• The probability of observing the pair of responses is determined by their frequencies in the population.
  • It’s the product of the frequences of those responses.
  • For ‘Yes,Yes’ it’s \(\theta_1^2\). For ‘No,No’ it’s \(\theta_0^2\).
  • And for ‘Yes,No’ and ‘No,Yes’ it’s \(\theta_1\theta_0\).
• This is, for what it’s worth, the product of the marginal probabilities of those responses.
  • When this happens for any pair of random variables, we say they are independent. More on that soon.
• And the probability of any response pair \(a,b\) depends only on its sum \(a+b\).

\[ \begin{aligned} \sum\limits_{\substack{j_1,j_2\\ y_{j_1},y_{j_2} = a,b}} \frac{1}{m^2} &= \sum\limits_{\substack{j_1 \\ y_{j_1}=a}} \qty{ \sum\limits_{\substack{j_2 \\ y_{j_2}=b}} \frac{1}{m^2} } \\ &= \sum\limits_{\substack{j_1 \\ y_{j_1}=a}} \qty{ m_b \times \frac{1}{m^2} } \qfor m_y = \sum\limits_{j:y_j=y} 1 \\ &= m_a \times m_b \times \frac{1}{m^2} = \theta_a \times \theta_b \qfor \theta_y = \frac{m_y}{m} \\ &= \theta_1^{s} \theta_0^{1-s} \qfor s = a+b \end{aligned} \]

• To calculate the marginal probability that \(Y_1+Y_2=s\), we sum the probabilities of all pairs with that sum.
  • This is pretty easy because there aren’t too many pairs.
  • And I’ve color coded the pairs to make it even easier.

• To find the distribution of a sum \(Y_1 + \ldots + Y_n\), we do the same thing.

• We start by writing out the joint distribution of \(n\) dice rolls.

  • We’re equally likely to roll any sequence of \(n\) numbers in \(1\ldots m\).
  • And there are \(m^n\) sequences, so the probability is \(1/m^n\) for each sequence.
  • Then we add columns for the corresponding responses ‘for free’. And the sums.

• Then we marginalize in two steps.

  1. Sum over roll sequences leading to the same response sequence \(a_1 \ldots a_n\).
  2. Sum over response sequences leading to the same sum \(s=a_1+\ldots+a_n\).

• This isn’t a class about counting, so this stuff won’t be on the exam.

• To sum over roll sequences, we sum over one roll after another.
  • This exactly like what we did for one response. Except over and over.
  • You can think of this sum ’inside out.
• Suppose we’re calculating the probability of the response sequence \(a_1 \ldots a_n = 0,0,...,0,0\).
  1. We consider any ‘head’ \(j_1\ldots j_{n-1}\) where \(y_{j_1} \ldots y_{j_{n-1}}=0,0,...0\) and think about choosing the last element. No matter what the head is, there are \(m_0\) ways to complete with by choosing \(j_{n}\) with \(y_{j_n}=0\).
  2. We repeat with shortened ‘head’ \(j_1 \ldots j_{n-2}\) to consider the choice of \(j_{n-1}\). We have \(m_0\) choices and our last step tells us that each choice results in \(m_0\) call sequences, so there are \(m_0^2\) sequences with head \(j_1 \ldots j_{n-2}\) ending in \(0,0\).
  3. Repeat \(n-2\) more times.
• What we get is pretty similar to the two-response case.
• The probability of observing the sequence is determined by the frequency of Yeses and Nos in the population.
  • It’s the product of the frequences of those responses.
  • For ‘Yes,Yes,…,Yes’ it’s \(\theta_1^n\). For ‘No,No,…,No’ it’s \(\theta_0^n\).
  • And if we have \(s\) Yeses and therefore \(n-s\) Nos, it’s \(\theta_1^s\theta_0^{n-s}\).
• This is, again, the product of the marginal probabilities of those \(n\) responses.

\[ \begin{aligned} \sum_{\substack{j_1 \ldots j_n \\ y_{j_1} \ldots y_{j_n} = a_1 \ldots a_n}} \frac{1}{m^n} &= \sum_{\substack{j_1 \ldots j_{n-1} \\ y_{j_1} \ldots y_{j_{n-1}} = a_1 \ldots a_{n-1}}} \sum_{\substack{j_n \\ y_n=a_n}} \frac{1}{m^n} \\ &= \sum_{\substack{j_1 \ldots j_{n-1} \\ y_{j_1} \ldots y_{j_{n-1}} = a_1 \ldots a_{n-1}}} \sum_{\substack{j_n \\ y_n=a_n}} m_{a_n} \times \frac{1}{m^n} \\ &= \sum_{\substack{j_1 \ldots j_{n-2} \\ y_{j_1} \ldots y_{j_{n-2}} = a_1 \ldots a_{n-2}}} \sum_{\substack{j_{n-1} \\ y_n=a_{n-1}}} m_{a_n} \times \frac{1}{m^n} \\ &= \sum_{\substack{j_1 \ldots j_{n-2} \\ y_{j_1} \ldots y_{j_{n-2}} = a_1 \ldots a_{n-2}}} \sum_{\substack{j_{n-1} \\ y_n=a_{n-1}}} m_{a_{n-1}} m_{a_n} \times \frac{1}{m^n} \\ &= m_{a_1} \ldots m_{a_{n}} \times \frac{1}{m^n} \qqtext{after repeating $n-2$ more times} \\ &= \theta_{a_i} \ldots \theta_{a_n} \\ &= \prod_{i:a_i=1} \theta_1 \prod_{i:a_i=0} \theta_0 = \theta_1^{s} \theta_0^{n-s} \qfor s = a_1 + \ldots + a_n \end{aligned} \]

• The full marginalization step is a bit trickier. But we can take advantage of our partial marginalization.
• To calculate the probability that the sum takes on the value \(s\), we order our sum of probabilities thoughtfully.
  • We sum the over response sequences \(a_1 \ldots a_n\) with \(s\) Yeses and \(n-s\) Nos.
  • And within that, sum over calls where we’d hear those responses.
• The inner sum we’ve done. That was our partial marginalization.
• And the outer sum boils down to counting the number of sequences with \(s\) Yeses.
  • The probability of all these sequences is the same: \(\theta_1^s\theta_0^{n-s}\). So summing is counting and multiplying.
  • You may have done the counting part in high school in a unit on ‘permutations and combinations’.
  • In any case, that count has a name. It’s spoken ‘\(n\) choose \(s\)’ and written \(\binom{n}{s}\).

\[ \begin{aligned} \sum\limits_{\substack{j_1 \ldots j_n \\ y_{j_1} + \ldots + y_{j_n} = s}} \frac{1}{m^n} &= \sum\limits_{\substack{a_1 \dots a_n \\ a_1 + \ldots + a_n = s}} \ \ \sum_{\substack{j_1 \ldots j_n \\ y_{j_1} \ldots y_{j_n} = a_1 \ldots a_n}} \frac{1}{m^n} \\ & =\sum\limits_{\substack{a_1 \dots a_n \\ a_1 + \ldots + a_n = s}} \theta_1^{s}\theta_0^{n-s} \\ &= \binom{s}{n} \theta_1^{s} \theta_0^{n-s} \qqtext{ where } \binom{s}{n} = \sum\limits_{\substack{a_1 \dots a_n \\ a_1 + \ldots + a_n = s}} 1 \end{aligned} \]

\[ P\qty(\sum_{i=1}^n Y_i = s) = \binom{s}{n} \theta_1^{s}\theta_0^{n-s} \ \ \text{ where } \ \ \binom{s}{n} \text{ is the number of binary sequences $a_1 \ldots a_n$ summing to $s$.} \]

• We don’t actually have to do all this sequence-summing-to-\(s\) counting ourselves.
• The choose function in R will do it for us: \(\binom{s}{n}\) is choose(n,s).
• The dbinom function will give us the whole probability: \(\binom{s}{n}\theta_1^s\theta_0^{n-s}\) is dbinom(s, n, theta_1).
• The rbinom function will draw samples from this distribution: rbinom(10000, n, theta_1) gives us 10,000.

theta_1 = .7
n = 625
p = dbinom(0:n, n, theta_1)
S = rbinom(10000, n, theta_1)

ggplot() + geom_area(aes(x=0:n,   y=p), color=pollc, alpha=.2) +
           geom_bar(aes(x=S, y=after_stat(prop)), alpha=.3) +
           geom_point(aes(x=S[1:100],   y=max(p)*seq(0,1,length.out=100)), 
                      color='purple', alpha=.2)

• To find the distribution of a sum \(Y_1 + \ldots + Y_n\) when we sample without replacement we do the same thing.

• We start by writing out the joint distribution of \(n\) dice rolls.

  • We’re equally likely to pull any sequence of \(n\) distinct numbers in \(1\ldots m\).
  • There are \(m \times (m-1) \times \ldots \times (m-n+1)=m!/(m-n)!\) sequences…
  • … so the probability is \((m-n)!/m!\) for each sequence.
  • Then we add the responses and sums as columns. That doesn’t change.

• Then we marginalize in two steps.

  1. Sum over roll sequences leading to the same response sequence \(a_1 \ldots a_n\).
  2. Sum over response sequences leading to the same sum \(s=a_1+\ldots+a_n\).

• To sum over roll sequences, we sum over one roll after another.
• Suppose we’re calculating the probability of the response sequence \(a_1 \ldots a_n = 0,0,...,0,0\).
  1. We consider any ‘head’ \(j_1\ldots j_{n-1}\) where \(y_{j_1} \ldots y_{j_{n-1}}=0,0,...0\) and think about choosing the last element. No matter what the head is, we’ve ‘used up’ \(n-1\) people with response \(y_j=0\). So there are \(m_0-(n-1)=m_0-n+1\) ways to complete the sequence by choosing an as-yet unused \(j_{n}\) with \(y_{j_n}=0\).
  2. We repeat with shortened ‘head’ \(j_1 \ldots j_{n-2}\) to consider the choice of \(j_{n-1}\). We have \(m_0-n+2\) as-yet unused choices and our last step tells us that each choice results in \(m_0-(n-1)\) call sequences, so there are \((m_0 - n + 1)(m_0 - n + 2)\) sequences with head \(j_1 \ldots j_{n-2}\) ending in \(0,0\).
  3. Repeat \(n-2\) more times. We get \((m_0 - n + 1) \times \ldots \times m_0 = m_0!/(m_0-n)!\) rows.
• When we work through the same argument for sequences with a mix of 0s and 1s, we see that a call in the head forecloses possibilities for later calls with the same response only. Thus, when we have \(s_1\) ones and \(s_0\) zeros, we can think of ourselves as choosing calls resulting with response zero and response one separately. We have \(m_0!/(m_0-s_0)! \times m_1! / (m_1-s_1)!\) rows.

\[ \begin{aligned} \sum_{\substack{j_1 \neq \ldots \neq j_n \\ y_{j_1} \ldots y_{j_n} = 0 \ldots 0}} \frac{m-n!}{m!} &= \sum_{\substack{j_1 \neq \ldots \neq j_{n-1} \\ y_{j_1} \ldots y_{j_{n-1}} = 0 \ldots 0}} \sum_{\substack{j_n \not\in j_1 \ldots j_{n-1} \\ y_n=0}} \frac{(m-n)!}{m!} \\ &= \sum_{\substack{j_1 \ldots j_{n-1} \\ y_{j_1} \ldots y_{j_{n-1}} = a_1 \ldots a_{n-1}}} \sum_{\substack{j_n \not\in j_1 \ldots j_{n-1} \\ y_n=0}} (m_{0} - n + 1) \times \frac{(m-n)!}{m!} \\ &= \sum_{\substack{j_1 \neq \ldots \neq j_{n-2} \\ y_{j_1} \ldots y_{j_{n-2}} = 0 \ldots 0}} \sum_{\substack{j_{n-1}\not\in j_1 \ldots j_{n-2} \\ y_n=0}} (m_0 - n + 2) (m_{0} - n + 1) \times \frac{(m-n)!}{m!} \\ &= m_0 \times \ldots \times (m_0 - n + 1) \times (m-n)!/m! = \frac{m_0!}{(m_0-n)!} \times \frac{(m-n)!}{m!} \qqtext{after repeating $n-2$ more times} \\ &= \frac{m_0!}{(m_0-s_0)!} \times \frac{m_1!}{(m_1-s)!} \times \frac{(m-n)!}{m!} \qqtext{ generally. } \end{aligned} \]

• The full marginalization step is exactly the same as it was for sampling with replacement.
  • Again, our response sequence probabilities depend only on the sequence’s sum.
  • So summing amounts to multiplying that probability by the number of sequences with a given sum.
  • And we counted those when we sampled with replacement. There are \(\binom{n}{s}\).

\[ \begin{aligned} \sum\limits_{\substack{j_1 \ldots j_n \\ y_{j_1} + \ldots + y_{j_n} = s}} (m-n)!/m! &= \sum\limits_{\substack{a_1 \dots a_n \\ a_1 + \ldots + a_n = s}} \ \ \sum_{\substack{j_1 \ldots j_n \\ y_{j_1} \ldots y_{j_n} = a_1 \ldots a_n}} (m-n)!/m! \\ & =\sum\limits_{\substack{a_1 \dots a_n \\ a_1 + \ldots + a_n = s}} \frac{m_0!}{(m_0-s_0)!} \times \frac{m_1!}{(m_1-s_1)!} \times \frac{(m-n)!}{m!} \\ &= \binom{s}{n} \frac{m_0!}{(m_0-s_0)!} \times \frac{m_1!}{(m_1-s)!} \times \frac{(m-n)!}{m!} \end{aligned} \]

\[ P\qty(\sum_{i=1}^n Y_i = s) = \binom{s}{n} \frac{m_0!}{(m_0-\{n+s\})!} \times \frac{m_1!}{(m_1-s)!} \times \frac{(m-n)!}{m!} \]

• This is called the hypergeometric distribution. R has what we need for this one too.
• The dhyper function will give us the probability of a sum \(s\). That’s dhyper(s,m_0,m_1,n).
• The rhyper function will draw samples: rhyper(10000, m_0, m_1, n) gives us 10,000.

n = 625
p = dhyper(0:n, sum(y==1), sum(y==0), n)
S = rhyper(10000, sum(y==1), sum(y==0), n)

ggplot() + geom_area(aes(x=0:n,   y=p), color=pollc, alpha=.2) +
           geom_bar(aes(x=S, y=after_stat(prop)), alpha=.3) +
           geom_point(aes(x=S[1:100],   y=max(p)*seq(0,1,length.out=100)), 
                      color='purple', alpha=.2)