Expectation of maximum of K Erlang variables: Am I doing it right?

Question

I want to find the expected value for the maximum of $K$ Erlang random variables (sum of two i.i.d exponential random variables with parameter $\gamma$), where $K = \frac{1}{\gamma}$. Below is what I have so far.

Let $$T = max(Z_1, Z_2, \cdots , Z_K)$$ where $Z_i$ is an Erlang random variable. The CDF of each Erlang random variable is $$F(X) = 1 - (1 + \gamma\cdot x)\cdot e^{-\gamma\cdot x}$$ The CDF for the maximum is then $$F_T(t) = P(T \leq t) = P(Z_1 \leq t, Z_2 \leq t, \cdots , Z_K \leq t)$$ Because all random variables are i.i.d, $$F_T(t) = P(Z_1 \leq t)\cdot P(Z_2 \leq t)\cdot~\cdots~\cdot P(Z_k \leq t) = F(t)^K = (1 - (1 + \gamma\cdot t)\cdot e^{-\gamma\cdot t})^K$$ To find the expected value, I must first differentiate $F_T(t)$ to find $f_T(t)$ which comes out to be $$f_T(t) = \gamma\cdot t\cdot e^{-\gamma\cdot t}\cdot[1 - (1 + \gamma t)\cdot e^{-\gamma\cdot t}]^{k - 1}$$ The expected value is then $$E[T] = \int_0^\infty t\cdot f_T(t)~dt = \gamma\int_0^\infty t^2\cdot e^{-\gamma\cdot t}\cdot[1 - (1 + \gamma t)\cdot e^{-\gamma\cdot t}]^{k - 1}~dt$$ Using the code below in MATLAB for $K = 20$ for numerical computation,

K = 20;
g = 1/K; % gamma
upper_bound = 10000
average = g * quad(@(x) f(x, K, g), 0, upper_bound)

function ret = f(x, k, g)
    ret = (x.^2).*exp(-g.*x).*((1 - exp(-g.*x).*(1 + g.*x)).^(k - 1));
end

I get

average =

 4.020119085379661e-08

which is definitely incorrect. If I increase upper_bound to say 1000000, the average drops to $0$. Using a different simulation for the same value of $K$, the average comes out to be 1.087141200000000e+02 which looks reasonable. What am I doing wrong in calculating the expected value, $E[T]$? The derivative of $F_T(t)$ seems to be fine.

Answer 1

It turns out that the problem was caused by the quad function. The documentation on this suggests to use the integral function instead. Once I switched quad with integral, I got the right answer. I am not sure if this is a bug.