Evaluation of the Laplace transform of the Pareto distribution

This function evaluates the Laplace-Stieltjes transform of the Pareto distribution

LAPLACEOFPARETO evaluates the Laplace-Stietljes transform (LST) of the Pareto
distribution with parameters a and b. Specifically, this function evaluates
the LST E[exp(-sX)] of a random variable $X$ which has Pareto distribution
$F(x)=1-\frac{b^a}{(b+x)^a}$, where s>0, a>0, b>0. The formula used can
be found in the following paper:
Nadarajah, S., Kotz, S. (2006) On the Laplace Transform of the Pareto
Distribution. Queueing Syst., 54, pp. 243-244.

SYNOPSIS: fun=LaplaceOfPareto(s,a,b)

INPUT s : argument
a,b : Pareto distribution parameters (positive reals)

OUTPUT fun : a vector of the same length as s; it contains NaN values
for zero values of the argument s. Obtain the value of the LST at zero
as a right limit, or develop this function further by using growth estimation of
the incomplete gamma function at zero.

REMARKS
1. This function makes use of the function gamma_incomplete:
http://www.mathworks.com/matlabcentral/fileexchange/29321-evaluation-of-the-upper-incomplete-gamma-function
The function gamma_incomplete.m evaluates the upper incomplete Gamma function and extends
the MATLAB function gammainc.
2. The Laplace-Stieltjes transform of a probability distribution is a
complete monotone function; it equals unity at s=0 iff the distribution
is proper (no positive mass at infinity), and it is less than unity at
s=0 if the distribution is improper; the difference 1-LST(0) is then the
mass at the infinity.


EXAMPLES

a=2.35; b=1.98;
h = 0.0001; s = h:h:100;
f = LaplaceOfPareto(s,a,b);
plot(s,f);

a=1; b=1;
h = 0.0001; s = h:h:200;
f = LaplaceOfPareto(s,a,b);
plot(s,f);

a=0.8; b=0.1;
h = 0.0001; s = h:h:500;
f = LaplaceOfPareto(s,a,b);
plot(s,f);

a=.1; b=0.1;
h = 0.0001; s = h/10:h:.1;
f = LaplaceOfPareto(s,a,b);
plot(s,f);

a=.1; b=0.1;
h = 0.0001; s = h/10:h:25;
f = LaplaceOfPareto(s,a,b);
plot(s,f);