What's the chance of drawing 8 or more cards of one suit if you draw 26 cards from a deck? My original question was going to be: "What's the chance of drawing 26 cards and not getting 7-7-6-6 cards of all four suits?" but I think the question "What's the chance of drawing 8 or more cards of one suit if you draw 26 cards from a deck?" is equivalent and might be more straightforward, I don't know.
I'd like to know how you'd go about solving a problem like this. I don't know how to start; my engineer friend told me that they'd use "formal methods" to solve this, but I can't wrap my head around it.
 A: The chance of getting $7,7,6,6$ is $4 \choose 2$ ways to choose the suits that get $7$ cards times ${13 \choose 7}^2$ ways to select those cards times ${13 \choose 6}^2$ ways to select the other cards out of ${26 \choose 13}$ total draws, so $$\frac {{4 \choose 2}{13 \choose 7}^2{13 \choose 6}^2}{26 \choose 13}\approx 0.105$$  Computing $7,7,7,5$ is similar.  Then to get the chance you don't get $8$ of a suit is the sum of these subtracted from $1$.
A: Let our sample space, $\Omega$, be all ways to draw twenty-six cards (where order doesn't matter).  We have then $|\Omega|=\binom{52}{26}$
Let $A_\spadesuit$ represent the event "drew at least 8 spades."  Similarly define $A_\clubsuit, A_\heartsuit, A_\diamondsuit$ as having drawn at least 8 clubs, hearts, and diamonds respectively.
We ask, what is $|A_\spadesuit\cup A_\clubsuit\cup A_\heartsuit\cup A_\diamondsuit|$?  I.e., how many ways can we have at least eight of some suit (possibly more than one).
Breaking apart via inclusion-exclusion, we have
$|A_\spadesuit\cup A_\clubsuit\cup A_\heartsuit\cup A_\diamondsuit| = |A_\spadesuit|+|A_\clubsuit|+|A_\heartsuit|+|A_\diamondsuit|-|A_\spadesuit\cap A_\clubsuit|-|A_\spadesuit\cap A_\heartsuit| - \dots + |A_\spadesuit\cap A\clubsuit\cap A\heartsuit| +\dots - |A_\spadesuit\cap\dots\cap A_\diamondsuit|$
Abusing the symmetry of the problem, we simplify as
$=4|A_\spadesuit| - 6|A_\spadesuit\cap A_\clubsuit| + 6|A_\spadesuit \cap A_\clubsuit\cap A_\heartsuit|$
noting that it is impossible to have at least 8 of every suit in hand at once (though it is possible to have at least eight of three suits at a time).
Calculating $|A_\spadesuit|$ is relatively straightforward as $\sum\limits_{k=8}^{13}\binom{13}{k}\binom{39}{26-k}$
$|A_\spadesuit\cap A_\clubsuit|$ is a bit messier.  $\sum\limits_{k=8}^{13}\sum\limits_{j=8}^{13}\binom{13}{k}\binom{13}{j}\binom{26}{26-k-j}$
Similarly, $|A_\spadesuit\cap A_\clubsuit\cap A_\heartsuit|$ will be $\sum\limits_{k=8}^{13}\sum\limits_{j=8}^{13}\sum\limits_{l=8}^{13}\binom{13}{k}\binom{13}{j}\binom{13}{l}\binom{26}{26-k-j-l}$
This final term could be trimmed a good deal, noting that if you have at least eight in three suits, you have either zero, one, or two cards of the off-suit, so the vast majority of terms in the above summation will be zero.
Regardless, using your favorite computing software, you could calculate the above number of cases and divide by the sample space size to get the probability.
(that being said, Ross's answer above is much easier to compute and turns out to be perhaps the more useful answer here due to the small number of cases in the complementary event.)
A: You'd use the hypergeometric distribution for drawing without replacement. We have $N = 52$ for the 52 total cards in the deck, $K = 13$ for the $13$ cards of a given suit, $n = 26$ for the number of draws, and $k \geq 8$ for the desired amount of successes. The probability mass function for a given number of successes is 
$$ 
\frac{\binom{K}{k}\binom{N - K}{n - k}}{\binom{N}{n}}
$$
To answer your question, we'd take the sum of this from $8$ to $13$ (because we want at least 8 successes, not exactly 8). Our expression for a single suit is
\begin{align*}
\sum^{13}_{k = 8}\frac{\binom{K}{k}\binom{N - K}{n - k}}{\binom{N}{n}} &= 
\sum^{13}_{k = 8}\frac{\binom{13}{k}\binom{52 - 13}{26 - k}}{\binom{52}{26}} \\
&= \sum^{13}_{k = 8}\frac{
\frac{13!}{(13 - k)!k!}
\frac{39!}{(13 + k)!(26 - k)!}
}
{
\frac{52!}{26!26!}
} \\
&= \sum^{13}_{k = 8} \frac {
13!39!26!26!
} {
52!(13 - k)!k!(k + 13)!(26 - k)!
} \\
&= \sum^{13}_{k = 8} \frac {26!26!} {52!}
\frac {13!} {k!(13 - k)!}
\frac {39!} {(26 - k)!(13 + k)!
} \\
&= \sum^{13}_{k = 8} \frac{1}{\binom{52}{26}}\binom{13}{k}\binom{39}{26 - k}
\end{align*} 
We then use the inclusion-exclusion principle, as detailed in JMoravitz's answer.
A: Ross Millikan's answer is a correct analytic answer, but it's nice to be able to work out the numerical answer via a simulation study in order to get more confidence in our answer.  In R this can be done in a few lines:
deck = rep(c("c","d","h","s"), 13)
set.seed(1)
table(replicate(10000, max(table(sample(deck, 26)))))

The variable deck represents just the suits that occur in the deck, so sample(deck, 26) is a sample of the 26 suits.  Applying table to that gives a table of the suits that occur, and the max of that table is the number of times the most frequent suit appears.  replicate replicates this 10000 times and table tabulates the results of that replication.
The return I get is
   7    8    9   10   11   12   13 
1552 4573 2823  890  146   15    1 

and so the largest suit has 7 cards (i. e. the 7-7-6-6 or 7-7-7-5 arrangement) 1552 times out of 10000 in my simulation; it has eight or more cards the rest of the time, 8448 times out of 10000.
