Is there a way to prove $\int {x^n e^x dx} = e^x \sum_{k = 0}^n {( - 1)^k \frac{{n!}}{{(n-k)!}}x^{n-k} } + C$ combinatorially? In How to integrate $\int x^n e^x dx$?, it is shown that
$$\int {x^n e^x dx}  = e^x \sum_{k = 0}^n ( - 1)^k \frac{n!}{(n-k)!}x^{n-k}   + C.$$
Since $\frac{n!}{(n-k)!}$ is $P(n,k)$, the number of $k$-permutations of $n$ objects, the right-hand side of this formula smells like the  inclusion-exclusion principle to me.  In fact, the right-hand side is almost the formula for the number of derangements on $n$ elements, which is one of the classic examples of the use of the inclusion-exclusion principle. 
This would seem to imply that there is some way to interpret $\int x^n e^x dx$ combinatorially such that the right-hand side formula falls out immediately.  However, I haven't been able to see it.  (I don't even know much about combinatorial interpretations of indefinite integration - or even if there are known ways to do so in special cases like this one.)  
Does anyone know or see a way to prove 
$$\int {x^n e^x dx}  = e^x \sum_{k = 0}^n ( - 1)^k \frac{n!}{(n-k)!}x^{n-k}   + C$$
combinatorially?
Added: To clarify, I am looking for a proof that shows that the left- and right-hand sides of the formula count the same thing.  That, to me, is a combinatorial proof.
 A: Depends what you mean by "combinatorics."  Does more analytic combinatorics and exponential generating series count?  Lets look at the definite integral $\int_0^x y^ne^y dy$ instead of $\int x^ne^xdx$
Let $a_n=\int_0^x y^n e^y$.  Then $$\sum_{n=0}^{\infty}a_{n}\frac{z^{n}}{n!}=\int_0^x e^{zy+y}dy=\frac{1}{\left(z+1\right)}e^{(z+1)x}-1=e^{x}\left(e^{zx}\left(1+z\right)^{-1}\right)-1.$$  
Taking the product of power series we have 
$$e^{x}\left(1+zx+\frac{\left(zx\right)^{2}}{2!}+\frac{(zx)^{3}}{3!}\cdots\right)\left(1-z+z^{2}-z^{3}+z^{4}+\cdots\right)-1.$$
$$=e^{x}\left(\sum_{n=1}^{\infty}\left(\sum_{k=0}^{n}(-1)^{n-k}\frac{x^{k}}{k!}\right)z^{n}\right)-1$$ and hence by comparing coefficients we conclude $$\int x^{n}e^{x}dx=n!e^x\sum_{k=0}^{n}(-1)^{n-k}\frac{x^{k}}{k!}+C.$$
Maybe that helps,
A: Though it is probably not what you are looking for, I find it relevant to give the following very simple derivation.
So we actually want to prove
$$
\frac{d}{{dx}}e^x \sum\limits_{k = 0}^n {( - 1)^k \frac{{n!}}{{(n - k)!}}x^{n - k} }  = x^n e^x .
$$
By the product rule, the left-hand side is given by
$$
e^x \sum\limits_{k = 0}^n {( - 1)^k \frac{{n!}}{{(n - k)!}}x^{n - k} }  + e^x \sum\limits_{k = 0}^{n - 1} {( - 1)^k \frac{{n!}}{{(n - k)!}}} (n - k)x^{n - k - 1} ,
$$
or
$$
e^x \bigg(x^n  + \sum\limits_{k = 1}^n {( - 1)^k \frac{{n!}}{{(n - k)!}}x^{n - k} }  + \sum\limits_{k = 0}^{n - 1} {( - 1)^k \frac{{n!}}{{(n - k - 1)!}}x^{n - k - 1} \bigg)} ,
$$
or
$$
e^x \bigg(x^n  + \sum\limits_{k = 1}^n {( - 1)^k \frac{{n!}}{{(n - k)!}}x^{n - k} }  + \sum\limits_{k = 1}^n {( - 1)^{k - 1} \frac{{n!}}{{(n - k)!}}x^{n - k} \bigg)} ,
$$
hence by
$$
e^x x^n.
$$
A: I managed to find a simple probability argument that's in the spirit of what I was looking for.
First, do a variable switch, with $t = -x$: $\int x^n e^x dx = (-1)^{n+1} \int t^n e^{-t} dt$.
Now, consider the definite version of the second integral, with an extra factor of $\frac{1}{n!}$:  $$\int_0^T \frac{t^n e^{-t}}{n!} dt.$$
This is the probability that the sum of $n+1$ independent exponential(1) random variables (which has an Erlang (or gamma) distribution) is less than $T$.  
Since interarrival times in a Poisson process are independent and exponentially distributed, this is also the probability that at least $n+1$ events have occurred by time $T$ in a Poisson process with rate 1.  If $N(T)$ denotes the number of events that occur in $[0,T]$ in a Poisson(1) process, we know that $P(N(T) = k) = T^k e^{-T}/k!$.  Thus 
$$\int_0^T \frac{t^n e^{-t}}{n!} dt = P(N(T) \geq n+1) = 1 - P(N(T) \leq n) = 1 - \sum_{k=0}^n P(N(T) = k) = 1 - \sum_{k=0}^n \frac{T^k}{k!} e^{-T}.$$
Switching back in terms of $x$ and multiplying through by $n!$, we have
$$\int x^n e^x dx = - (-1)^{n+1} \sum_{k=0}^n \frac{n!}{k!} (-x)^ke^x +C = e^x \sum_{k=0}^n (-1)^{n-k} \frac{n!}{k!}x^k +C$$
$$= e^x \sum_{k=0}^n (-1)^k \frac{n!}{(n-k)!}x^{n-k} + C.$$
A: Here is a solution I quite like which uses inclusion exclusion on a integral over an $n$-dimensional simplex.
First, it is much nicer to consider the definite integral $\int_0^x y^n e^{y} dy$ which will only be off by a constant.  (Which we don't care about anyway)
Consider
$$\int_{0}^{x}\int_{0}^{t_{1}}\cdots\int_{0}^{t_{n}}e^{-t_{n+1}}dt_{n+1}dt_{n}\cdots dt_{1}.$$
This is $e^{-t_{n+1}}$ integrated over the $n$-dimesional simplex. (Only one coordinate is integrated) By Cauchy's formula for repeated integration, this is the same as $$\frac{1}{n!}\int_{0}^{x}(x-y)^{n}e^{-y}dy=\frac{1}{n!}e^{-x}\int_{0}^{x}y{}^{n}e^{y}dy.$$ Now, we need only show that $$\int_{0}^{x}\int_{0}^{t_{1}}\cdots\int_{0}^{t_{n}}e^{x-t_{n+1}}dt_{n+1}dt_{n}\cdots dt_{1}=e^{x}\sum_{k=0}^{n}\frac{(-1)^{n-k}}{k!}x^{k}+C.$$  (A constant will arise since the integrand is not $0$ at $0$) But this follows from splitting up the simplex and using inclusion exclusion.  Here is why:  Integrating $e^{x-t_1}$ over the cube gives $x^ke^x$.  Since the cube can be divided into $k!$ simplicies we get $\frac{x^k}{k!}e^x$.  Now we need to subtract off a correction factor based on the $n-1$-dimensional simplicies, and so on.
Hope that helps,
