Help with the integral for the variance of the sample median of Laplace r.v. When we draw $n$ samples of Laplace-distributed random variable such that $n=2k+1$ and the location parameter is zero, the median $x$ (or the $k$-th order statistic) has the following p.d.f.:
$$f_m(x)=\frac{n!}{(k!)^2}2^{-n}\frac{1}{b}e^{-(k+1)|x|/b}(2-e^{-|x|/b})^k$$
where the p.d.f. of the underlying Laplace distribution is given as $f(y)=\frac{1}{2b}e^{-|y|/b}$.
The formula for p.d.f. of the median stems from the usual method of characterizing the distributions of order statistics and is found as equation (2.5.10) in Kotz's volume on Laplace distribution.  There is another formula for the case when $n$ is even, but we shall not be concerned with it for now.
I am interested in the variance of the sample median.  Since $f_m(x)$ is symmetric about $x=0$, I can get rid of the absolute value and write it as follows:
$$\sigma^2_m=\frac{n!}{(k!)^2}2^{-n}\frac{1}{b}\int_{0}^{\infty}x^2e^{-(k+1)x/b}(2-e^{-x/b})^kdx$$
What I need is to reduce the following definite integral to a more manageable form:
$$\int_{0}^{\infty}x^2e^{-(k+1)x/b}(2-e^{-x/b})^kdx$$
For fixed $k$ this is fairly easy integral to solve.  However, I am interested in asymptotics of variance $\sigma_m^2$ as $n\rightarrow\infty$ and $b$ is a linear function of $\sqrt{n}$, and thus need a solution for an arbitrary $k$.
Any ideas?
 A: Ok, I think I understand now. You want to find
$$\sigma^2=\frac{n!}{(k!)^2}2^{-n}\frac{1}{b}\int_{0}^{\infty}\!dx\,x^2e^{-(k+1)x/b}(2-e^{-x/b})^k$$
where $n=2k+1$ and $b \sim \sqrt{n}$ for $n,k \to \infty$. Let us first use the change of variable $y=x/b$. We then find
$$\sigma^2= b^2\frac{n!}{(k!)^2}2^{-n}\int_{0}^{\infty}\!dy\,y^2e^{-(k+1)y}(2-e^{-y})^k.$$
We will use the method of steepest decent to figure out the asymtptotic behavior of $\sigma^2$ for $k,n \to \infty.$ Let us first rewrite the integral as
$$ \int_0^\infty\!dy\,y^2e^{-(k+1)y}(2-e^{-y})^k =
\int_0^\infty\!dy\, y^2 e^{-y} e^{- k f(y)}, \qquad f(y) =y-\log(2-e^{-y}). $$
The function $f(y)$ assumes its maximum at $y=0$. Thus, the integral is asymptotically dominated for values of $y$ close to $0$. To this end, we expand $f(y) = y^2 +O(y^3)$ and we  have
$$\int_0^\infty\!dy\, y^2 e^{-y} e^{- k f(y)} \sim
\int_0^\infty\!dy\, (y^2 -y^3 +O(y^4)) e^{-k y^2} = \frac{\sqrt\pi}{4 k^{3/2}} - \frac{1}{2k^2} +O(k^{-5/2}). $$
Additionally, we use that (see central binomial coefficient) $$\frac{n!}{(k!)^2} \sim2 \binom{2k}{k} \sim 2 \frac{4^k}{\sqrt{\pi k}} = \frac{2^n}{\sqrt{\pi k}}.$$
Together, we have (with $b \sim \alpha \sqrt{n} \sim 2 \alpha \sqrt{k} $)
$$\sigma^2 \sim \underbrace{4 \alpha^2 k}_{b^2} \frac{2^n}{\sqrt{\pi k}} 2^{-n} \frac{\sqrt\pi}{4 k^{3/2}} = \frac{\alpha^2}{k}. $$
A: Let's first set $b=1$, and recover it later using scaling. For $X \sim \mathrm{Lap}(0,1)$
$$
   f_{X_{k:2k+1}}(x) = \frac{k+1}{2^{2k+1}} \binom{2k+1}{k+1} \mathrm{e}^{-(k+1)|x|} \left(2-\mathrm{e}^{-|x|}\right)^k = \frac{k+1}{2^{2k+1}} \binom{2k+1}{k+1} \sum_{m=0}^k \binom{k}{m} 2^{k-m} (-1)^m \exp\left(-(m+k+1)|x|\right)
$$
Thus
$$
  \mathbb{Var}\left(X_{k:2k+1}\right) = \frac{k+1}{2^{2k+1}} \binom{2k+1}{k+1} \sum_{m=0}^k \binom{k}{m} 2^{k-m} (-1)^m \frac{4}{(k+1+m)^3} = \\ \binom{2k+1}{k+1} \frac{2^{1-k}}{(k+1)^2} \cdot {}_4F_3\left( \left.\begin{array}{cccc} -k & k+1 & k+1 & k+1 \\
   & k+2 & k+2 & k+2 \end{array} \right| \frac{1}{2} \right)
$$

The sequence above satisfies an inhomogeneous recurrence equation of rank $3$:

This recurrence equation allows to determine the large $k$ asymptotic behavior. There are fewer undetermined coefficients, because $X_{k:2k+1}$ converges to a degenerate random variable for large $k$, meaning that the variance must vanish in this limit:
$$
   \mathbb{Var}\left(X_{k:2k+1}\right) = \frac{c_1}{k}\left(1+ \mathcal{o}\left(k^{-1}\right) \right) + \frac{c_2}{k^{3/2}}\left(1+ \mathcal{o}\left(k^{-1}\right) \right) 
$$
The following plot indicates that $c_1>0$:

Sequence acceleration  suggests that value of $c_1$ is close to $\frac{1}{2}$:

