Are the error terms of the partial sums of inclusion-exclusion unimodal?

The error terms need not be unimodal. Let $x$ be a point which is in $A_i$ for exactly $n \geq 1$ indices $i \in [k] := \{1,\ldots,k\}$. Knowing $n$, it's not actually that hard to give a formula for $x$'s contribution to the error term $R_m$ and I think that this makes what's going on pretty clear. Actually I will use $R_m$ for

$$ \left| \bigcup_{i=1}^k A_i \right| - \left( \sum_{r=1}^m (-1)^{r+1} \sum_{I \in \binom{[k]}{r}} \left| \bigcap_{i \in I} A_i \right| \right)$$

ie. the signed error which results when we count only the contribution to the inclusion-exclusion formula from $r$-fold intersections, $1 \leq r \leq m$.

Since $n \geq 1$, $x$ should contribute exactly $1$ to $| \bigcup_{i =1}^k A_i |$. For positive integers $r$, the number of $r$-element subsets $I$ of $[k]$ for which $x \in \bigcap_{i \in I} A_i$ is $\binom{n}{r}$ (choosing amoung the $n$ indices $i$ for which $x \in A_i$). Thus $x$'s contribution to $R_m$ is simply

$$ 1 - \left[ \binom{n}{1} - \binom{n}{2} + \ldots \pm \binom{n}{m} \right] = \sum_{r = 0}^m (-1)^r \binom{n}{r} = (-1)^m \binom{n-1}{m}$$

where the closed form can be obtained using the identity $\binom{a+1}{b+1} = \binom{a}{b} + \binom{a}{b+1}$ and induction. Unfortunately I can't see a combinatorial proof. Now it is clear that $x$'s contribution to $R_m$ is most significant when $m \sim (n-1)/2$ and is $0$ for $m \geq n$.

It's easy to disprove unimodality now. For example, if we set things up so that the universe consists of a large number $N$ of points which are in $A_i$ for say $10$ of the $i \in [k]$ and a single further point which is in $A_i$ for $100$ of the $i \in [k]$, then $|R_m|$ will have one a peak around $m=4,5$ when the $N$ points have the most effect. The additional point will contribute the most to the error when $m \sim 50$ and this will result in another peak in $R_m$ since the $N$ points will have long since stopped having any effect after $m=10$.

Added: As Gerry's answer shows, my proposed counterexample above is quite suboptimal. I hadn't played around with the numerics at all when I gave it so I made the numbers big out of some vague paranoia about edge effects.

I agree that Gerry's example simultaneously minimizes the number of sets in the collection and the number of elements in the universe. If $c_n \geq 0$ is the number of elements in exactly $n$ sets of the collection for $n=1,2,\ldots,k$ then, from the argument above, we have

$$R_m = (-1)^m \sum_{n=1}^k c_n \binom{n-1}{m}$$

So, the possibilities for the sequence of absolute remainders are exactly the nonnegative integer linear combinations of the rows of the array:

$$\begin{array}{clcr} 0 & 0 & 0 & 0 & 0 & 0 & 0 & \cdots \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & \cdots \\ 2 & 1 & 0 & 0 & 0 & 0 & 0 & \cdots \\ 3 & 3 & 1 & 0 & 0 & 0 & 0 & \cdots \\ 4 & 6 & 4 & 1 & 0 & 0 & 0 & \cdots \\ 5 & 10 & 10 & 5 & 1 & 0 & 0 & \cdots \\ 6 & 15 & 20 & 15 & 6 & 1 & 0 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \end{array}$$

For example, the absolute remainders from Gerry's example are obtained by taking nonzero coefficents 6, 4 and 1 for the 2nd, 3rd and 7th rows. I think that viewed this way it is perhaps easier to convince oneself that this is the "minimal" choice of coefficients such that the resulting linear combination exhibits non unimodality.

Here's an explicit example of non-unimodality. We use $7$ sets: $$\eqalign{A_1&=A_2=\lbrace a,b,c,d,e,f,g,h,i,j,k\rbrace\cr A_3&=\lbrace g,h,i,j,k\rbrace\cr A_4&=A_5=A_6=A_7=\lbrace k\rbrace\cr}$$ Then the absolute remainders go $20,19,20,15,6,1,0$.

This may be, in some sense, the minimal example. I don't think it can be done with fewer than $7$ sets, or fewer than $11$ elements all told (but my reasons involve considerable hand-waving and appeals to plausibility).

This isn't directly an answer to your question on unimodality, but I hope it's ok to post it since it's a pretty good answer to the question implicitly raised in the first part of your post - to what extent is it the case that final terms in the formula are "corrections" to the partial sum?

This paper of Linial and Nisan discusses the behavior of the sequence $R_m$ in terms of how closely it approximates the final answer. It turns out that the approximation is poor through $m=O(\sqrt{k})$, and then it starts to provide pretty good estimates of the final answer. I'm quite sure that the behavior in the first phase is not necessarily monotonic, meaning the answer to your question is "no", but I don't think this is directly addressed in the paper.