Conceptual understanding of waiting time problem involving exponential distributions I am working on a textbook problem that says,

A doorman at a hotel is trying to get three taxicabs for three different couples. The arrival of empty cabs has an exponential distribution with mean 2 minutes. Assuming independence, what is the probability that the doorman will get all three couples taken care of within 6 minutes?

(Robert V. Hogg et al., Probability and Statistical Inference 9e, exercise 5.4.19.)
I can see two ways of approaching this problem:
Each taxi's arrival time has the cdf $F(x) = 1-e^{-x/2}$. Then the probability of one given taxi arriving in 6 minutes is $1- e^{-6/2}$, and the probability of this happening three times is $\left(1- e^{-3}\right)^3 = 0.8580$ (answer 1).  This is equivalent to defining $Y = \max(X_1,X_2,X_3)$ and finding $P(Y \leq 6)$. 
But we could also say that $f(x)$ is counting the time until the "next" taxi. Then we define $X_1 + X_2 + X_3 = Z$, the total time until three taxis have shown up. Since each $X_i$ has the mgf $M_X(t) = (1-2t)^{-1}$,
$$M_Z(t) = \prod_{i=0}^3 M_{X_i} = (1-2t)^{-3}$$
This is gamma with $\alpha = 3$ and $\theta = 2$. Using Python, I get this (answer 2):
In [265]: from scipy.stats import gamma
In [266]: gamma.cdf(6, 3, 2)
Out[266]: 0.7618966944464556

The solutions page says, without explanation, that the answer is $ 1 −( 17/2 )e^3 = 0.5768$ [sic: should have $e^{-3}$] (answer 3).
What are the assumptions behind these three different approaches? Where did the textbook's answer come from? Which approach is correct?
 A: First, let us address the textbook's solution and its origin.  This solution relies on observing that the taxicab arrivals constitute a Poisson process with intensity $\lambda = 1/2$ taxis per minute.  Consequently, we seek the probability that by time $t = 6$ minutes, we have observed at least $3$ taxis.  The counting random variable for the number of taxis observed by time $t$ is $$X_t \sim \operatorname{Poisson}(\lambda t),$$ with $$\Pr[X_t = x] = e^{-\lambda t} \frac{(\lambda t)^x}{x!}, \quad x = 0, 1, 2, \ldots.$$  Consequently $$\Pr[X_6 \ge 3] = 1 - \Pr[X_6 \le 2] = 1 - e^{-3} \frac{3^0}{0!} - e^{-3} \frac{3^1}{1!} - e^{-3} \frac{3^2}{2!} = 1 - \frac{17}{2}e^{-3},$$ as claimed.
Now we turn our attention to the first solution you describe.  The flaw is one of interpretation.  The arrival of taxicabs is a serial process, not a parallel one; thus, the statement that the time of arrival of empty cabs is exponentially distributed means that these are interarrival times, not absolute times.  In other words, conditioned on the random arrival of a cab, the arrival time of the next cab is exponential with mean $2$ minutes.  The model used in the first solution's calculation might apply to a scenario where, for instance, three cabs are known to be arriving independently of each other, each one's arrival time is exponential with mean $2$ minutes, and you want to know the probability that the last cab to arrive will do so within $6$ minutes.  But this is not how the question is set up.  In particular, this solution fails to model situations where four cabs, or even ten cabs, might show up in $6$ minutes.
For the second solution, the method is correct, and you should have gotten the same numeric result as with the third (textbook) solution.  The time-to-event random variable $T_x$ that measures the time it takes to observe the $x^{\rm th}$ arrival is gamma with shape $x$ and intensity (rate) $\lambda = 1/2$, i.e. $$T_3 \sim \operatorname{Gamma}(3, \lambda = 1/2), \quad f_{T_x}(t) = \frac{\lambda^x t^{x-1} e^{-\lambda t}}{\Gamma(x)}, \\
\Pr[T_3 \le 6] = \int_{t=0}^6 \frac{t^2 e^{-t/2}}{16} \, dt.$$  I suspect you obtained the result you did because you specified the arguments incorrectly; however, I am not familiar with Python.  I believe you must type instead
gamma.cdf(6, 3, scale=2)

because otherwise Python will interpret the third argument as a location parameter, and leave the scale parameter at the default value of $1$.  So your syntax is actually equivalent to
gamma.cdf(6, 3, loc=2, scale=1)

which in turn is also equivalent to
gamma.cdf(4, 3) 

