I want to use this:
Let $A_1, A_2, \ldots, A_n$ be $n$ sets in a universe $U$ of $N$ elements. Let $S_1 = \sum_jN(A_j), \ S_2 = \sum_{i, j}N(A_iA_j), \ldots, S_k = \sum N(A_{j_1}A_{j_2\ldots A_{j_k}}), \ldots, S_n = N(A_1A_2\ldots A_n)$
Then $$ N(\overline A_1\overline A_2\ldots \overline A_n) = N(U) - S_1 + S_2 - S_3 + \ldots + (-1)^kS_k+ \ldots + (-1)^nS_n$$
to derive this:
The number of onto functions $\{1, 2, 3, \ldots, r\} \to \{1, 2, 3, \ldots, n\}$ is $$\sum_{j=0}^n\binom nj(-1)^j(n - j)^r$$.
So,
Let $S_0 = N(U).$ Then $$N(U) - S_1 + S_2 - S_3 + \ldots + (-1)^kS_k+ \ldots + (-1)^nS_n \\ = \sum_{k=0}^n(-1)^kS_k \\ = \sum_{k=0}^n(-1)^k\sum N(A_{j_1}A_{j_2\ldots A_{j_k}}) \\ = \sum_{k=0}^n(-1)^k \binom nkN(A_{j_1}A_{j_2\ldots A_{j_k}}) \\ = \sum_{k=0}^n(-1)^k \binom nkN(A_{j_k})$$
Now let $A_{j_k}$ be the set of functions that miss $k$ elements. Then $N(A_{j_k})$ is the number of such functions which is given by $(n - k)^r.$
The derivation above seems to work iff $N(A_{j_1}A_{j_2\ldots A_{j_k}}) = N(A_{j_k})$ which I am not sure is true. If so, how do we prove that? Thanks.
Notice that $N(A_{j_k})$ has $(n-1)^r$ because you are missing just $1$ element. Instead, $N(A_{j_1}A_{j_2}\cdots A_{j_k})$ is $(n-k)^r$ because you are missing now $k$ elements instead. They are not equal is just that the last equality is not true.