Is left-of-mean cdf of exponential distribution greater than that of a gamma distribution? Edited to restrict one of the gamma distributions to be an exponential distribution, and to force both means to 1.0, and to add wolfram alpha links:
I have two distributions, both with mean 1, and I want to compare their cdfs when $0 \lt x \lt 1$.  One of the distributions is exponential, and the other distribution is gamma.  I think that in the $0 \lt x \lt 1$ domain, the cdf of the exponential distribution is always at least as big as that of the gamma distribution, but I have not yet been able to prove it.
The exponential distribution of interest is completely specified by its mean of 1.0, and it can be seen for example here or here using wolfram alpha.
The gamma distribution of interest is forced to have mean 1.0 but it has an extra degree of freedom compared to the exponential distribution.  I want the gamma distribution to be the distribution with mean 1.0 that you get by summing $k \ge 1$ i.i.d exponential distributions each with mean $1 / k$.  Using wolfram alpha's notation this can be seen here (broken link).  Notice that when $k=1$ the gamma distribution is equal to the exponential distribution with mean 1.0, and that in general its cdf is $\Pr(X \le x) = Q(k, 0, k x)$.
This wolfram alpha plot (broken link) shows that for a few $k$ and for $x \in (0,1)$, the cdf seems to increase as a function of $x$ and to decrease as a function of $k$.  I would like to be able to prove the monotonicity in $k$; the monotonicity in $x$ is just a result of the definition of a cdf.  Note that the monotonicity in $k$ does not continue (broken link) when $x \ge 1$.
 A: The answer to your monotonicity conjecture is affirmative. It admits a somewhat sneaky proof and allows us to conclude something about the Poisson distribution in the
process. This is what we explore below.
The picture
The question asks whether the cdfs decrease pointwise as $n$ increases
for each $x$ to the left of the vertical line $x=1$.

The math
Let's start with a restatement of the question. First, if $\renewcommand{\Pr}{\mathbb P}\newcommand{\d}{\mathrm d}T \sim
\mathrm{Exp}(1)$, then it's easy to see that $T/n \sim
\mathrm{Exp}(n)$. So, we can rewrite a $\Gamma(n,n)$ random variable
$S_n$ as $S_n = (T_1 + \cdots + T_n)/n$ where $T_i$ is a sequence of
iid $\mathrm{Exp}(1)$ random variables. Instead of asking whether the
cdf decreases pointwise as a function of $n$ for each $0 < x < 1$, we
can ask whether the survival distribution increases
pointwise. Hence, we have the following equivalent question.

Equivalent question: Let $T_1, T_2,\ldots$ be an iid sequence of   $\mathrm{Exp}(1)$ random variables and let $S_n = T_1 + \dots +  T_n$. Then, is it true that, for all fixed $x \in (0,1)$, $$ \Pr( S_{n+1} > (n+1) x ) \geq \Pr(S_n > n x) \,? $$

Proof. The key idea here is to divide and conquer. We need to somehow compare the events $\{S_{n+1} > (n+1)x\}$ and $\{S_n > nx\}$. To that end, notice that
$$
\Pr(S_{n+1} > (n+1)x) = \Pr(S_{n+1} > n x) - \Pr(nx \leq S_{n+1} \leq (n+1)x) \>.
$$
We are halfway there. To deal with the remaining part, observe that
$$
\{S_{n+1} > n x\} = \{S_n > n x\} \cup \{S_n \leq n x, T_{n+1} > nx - S_n \} \>,
$$
and the two events on the right-hand size are disjoint. So,
$$
\Pr(S_{n+1} > (n+1)x) = \Pr(S_n > n x) + \Pr(S_n \leq nx, T_{n+1} > nx - S_n) - \Pr(nx \leq S_{n+1} \leq (n+1)x) \>.
$$
If we can show the second probability is greater than the third, we are done. To do this, we'll first find an explicit value for the second probability and then find an upper bound for the third.
(Second term.) $S_n$ and $T_{n+1}$ are independent, so to find $\Pr(S_n \leq nx, T_{n+1} > nx - S_n)$, we need only integrate over the
region of interest. This is
$$
\Pr(S_n \leq nx, T_{n+1} > nx - S_n) = \int_0^{nx} \int_{nx-s}^\infty
e^{-t} \frac{s^{n-1}e^{-s}}{(n-1)!} \, \d t \, \d s  =
\frac{e^{-nx}(nx)^n}{n!} \>.
$$
(Third term.) This step is more involved, but still employs only elementary tools. Using the density of $S_{n+1}$, we have
$$\newcommand{\d}{\mathrm d}
\Pr(nx \leq S_{n+1} \leq (n+1)x )= \int_{nx}^{(n+1)x} \frac{u^n
e^{-u}}{n!} \,d u = \frac{e^{-nx} (nx)^n}{n!} x \int_0^1 (1+v/n)^n
e^{-xv} \d v
$$
where we've obtained the integral on the far right by using the substitution $v =
(u-nx)/x$. Now $(1+v/n)^n < e^v$ for all $n$ and so
$$
\int_0^1 (1+v/n)^n e^{-xv} \d v < \int_0^1 e^{(1-x) v} \d v =
\frac{e^{1-x} - 1}{1-x} \leq \frac{1}{x} \>,
$$
where the last step follows from the fact that $e^{1-x} < x^{-1}$ for
all $x \in (0,1)$. So,
$$
\Pr(nx \leq S_{n+1} \leq (n+1)x ) < \frac{e^{-nx}(nx)^n}{n!} \>. 
$$
But, this last quantity is the same as the second term! Putting the three pieces together, we get
$$
\Pr(S_{n+1} > (n+1)x ) > \Pr(S_n > n x) \>,
$$
for each $x \in (0,1)$, which is what we set out to prove.
A consequence
The appearance of the term $\frac{e^{-nx}(nx)^n}{n!}$ strongly hints
that there is a relation to the Poisson hiding in here. Indeed, the
proof implies the following.

Proposition: Let $X_{n\lambda} \sim \mathrm{Poisson}(n\lambda)$   for $\lambda \in (0,1)$. Then $$ \Pr(X_{n\lambda} < n) \geq e^{-\lambda} \>. $$

Indeed, this was the original context in which I recently gave an equivalent
proof in this MathOverflow answer.
A: An analytical proof might be quite difficult to obtain. The CDF of a gamma distribution fixed to have mean 1 is $P(k,kx) = \frac{\gamma\left(k,kx\right)}{\Gamma\left(k\right)}$ and you want to show for $0 \le x \le 1$ and $k > 0$ that
$$
\frac{dP\left(k,kx\right)}{dk} \ge 0
$$
But this works out as being fairly horrible:
$$
\frac{dP\left(k,kx\right)}{dk} = - \frac{x \operatorname{\gamma}\left(k, k x\right) \operatorname{\psi}\left(k x\right)}{\operatorname{\Gamma}\left(k x\right)} + \frac{\frac{x \left(k x\right)^{k -1}}{e^{k x}} - \operatorname{log}\left(k x\right) \operatorname{\Gamma}\left(k, k x\right) + \operatorname{\Gamma}\left(k\right) \operatorname{\psi}\left(k\right) + {G_{2, 3}^{3, 0}\left.\left(\begin{matrix} 1, 1 \\0, 0, k  \end{matrix} \right| {k x} \right)}}{\operatorname{\Gamma}\left(k x\right)}
$$
where $\psi(x)$ is the digamma function and $G_{p,q}^{m,n}\left.\left(\begin{matrix} a_1,\ldots,a_p \\ b_1,\ldots, b_q \end{matrix} \right| z \right)$ is the Meijer G-function.
I'm not saying it's impossible to prove the inequality but it looks like an awful lot of work!
FYI the differential was calculated using Python 2.7 and the SymPy library.
PS there may be other ways to prove this, for example you could try expanding $P(k,kx)$ into its integral form and working from there
