Lambart function

[SWEETSAS]

Hello all:

Does anyone know whether any SAS procedure has Lambart function?

Thanks

J

[Rick_SAS]

I assume you mean the real-valued Lambert W function? No, I don't think that is built into Base SAS or SAS/IML.

To understand your needs better, could you share

1) how you hope to use the Lambert W function? What application do you have in mind?

2) The Lambert W function isn't actually a function. It has two branches. Do you need both branches?

Although the function isn't built in, you can implement it to arbitrary precision by using the fact that the function is the inverse of f(w)=w*exp(w).  It is easy to evaluate f.  For example, the following statements plot the upper branch of the Lambert W function:

proc iml;
w = do(-1, 2, 0.1);
x = w#exp(w);
title "Lambert's W Function";
call series(x, w) grid={x y};   /* inverse of w*exp(w) */

Since the function is monotonic on each branch, you can use root-finding to evaluate the inverse. 

[SWEETSAS]

Thanks Rick. 

Yes. I want to it use for generating random variates from a new function with a specified parameter.

I am generating data from a function that I created,

E.g: f(x)=theta*exp(-theta*x). I want to generate random variable from f(x). When theta is given say 0.1. 

I like to specifically use Lambert function for the generation

Thanks,

J

[Rick_SAS]

OK, so let's talk about generating random variates. 

The density function f(x) = theta*exp(-theta*x)

is an exponential distribution with "rate parameter" theta.

An alternate parameterization is exp(-x/sigma) / sigma, where sigma=1/theta is called the shape parameter.

The RANDGEN function in SAS/IML supports random generation from the exponential distribution, although it uses the shape parameter.  Consequently, here is how you can generate 500 random variates from the exponential distribution with rate parameter 0.1 in SAS/IML:

proc iml;
call randseed(12345);  /* set random number seed */
x = j(500, 1);         /* allocate a vector to hold results */
theta = 0.1;           /* rate parameter */
call randgen(x, "Exponential", 1/theta); /* shape parameter */

call histogram(x) rebin={2.5 5}; /* visualize the random sample */

[SWEETSAS]

Thank you very Rick.

I use SAS for my random number generation.

Actually, I am trying to convince some colleagues who is are fan of another software. We are developing a new density that has not been listed. I am recommending we use SAS. They have been using Lambert function implemented in another software for genetation since the density is not supported by software. Since Lambert function is the method they have used for many years for generating data from unlisted, I want to demonstrate that it can be done in SAS.

And show that the variate generatedd using Lambert function is no differrent that generated using inverse CDF solve using numerical methods.

Thanks

J

[Rick_SAS]

Sounds great.

• I have written about how to use the inverse CDF method to generate random variates.
• You can use the FROOT function in SAS/IML to numerically solve inverse problems.
• Chapter 7 of Simulating Data with SAS contains other advanced methods for simulating data from univariate distributions. Including root-finding as part of the invrse CDF method.

If you write out the PDF or CDF that you want, I might be able to suggest some efficient computational techniques.

[SWEETSAS]

Here are the CDF and the Qauntile function expressed in terms of Lambert function.

The goal is to generate the random variables t using1)  inverse CDF and 2) Lambert function method.

Distribution function: 

alphat=(alpha*t/(1+alpha));

expalpha=exp(-alpha*t);

F(t|alpha)=1-(1+alphat)*expalpha;

Quantile function (Q):

Q(P|alpha)=-1-(1/alpha)-(1/alpha)W(-(1+alpha)*(1-p)*exp(-(1+alpha));

0<p<1; Where W is the Lambert function. P is essentially a unif(0,1)

I will also like to track the computer time using both method. I want SAS to win

Thanks,

J

[Rick_SAS]

If your goal is to generate random variates for this distribution, I don't think you need to implement the inverse CDF method or use the Lambert W function. My calculations indicate that your density is a mixture distribution of two well-known densities: an exponential distribution and a gamma distribution. (Actually Erlang, which is a gamma distribution with an integer shape parameter.)

Here are my calculations: Instead of the CDF, look at the PDF. Then the density function is

f(t) = alpha^2 / (1+alpha) exp(-alpha*t)(1+t)

     = alpha/(1+alpha) [ alpha*exp(-alpha*t) ] +  1/(1+alpha) [ alpha^2 *t*exp(-alpha*t) ]

     = alpha/(1+alpha) [ g(t) ] +  1/(1+alpha) [ h(t) ]

where g(t) is the density of the exponential distribution with shape parameter sigma=1/alpha, and 

where h(t) is the density of the gamma distribution with integer shape parameter a=2 and second shape parameter sigma=1/alpha 

The following SAS/IML program compares your density to the mixture density:

proc iml;
alpha = 1.5;
t = do(0, 7, 0.1);
f = alpha##2 /(1+alpha) # (1+t) # exp(-alpha#t); /* mystery distrib */

shape = 1/alpha;
g = pdf("Expon", t, shape);    /* exponential distrib */
h = pdf("Gamma", t, 2, shape); /* gamma distrib (actually Erlang with a=2 */
p = alpha / (1+alpha);         /* mixing probability */
mixDist = p*g + (1-p)*h;

/* show that mixDist = mystery distrib */
maxDiff = max(abs(f - mixDist));
print maxDiff;

The densities are identical to machine precision.

By linearity, the PDF and CDF for your distribution is just a linear combination of the PDF/CDF for the exponential and gamma distributions. This blog post shows how to generate a random sample from a mixture distribution:

/* random sample from mixture distribution */
/* See http://blogs.sas.com/content/iml/2011/09/21/generate-a-random-sample-from-a-mixture-distribution.html */
N = 1000;
call randseed(12345);
Type = 1 + randfun(N, "Bernoulli", p); /* fill with 1s and 2s */

x = j(N,1); /* allocate vector */
/* Exponential draw */
idx = loc(Type=1);
y = j(ncol(idx),1);
call randgen(y, "Expo", shape);
x[idx] = y; 
/* Erlang (gamma) draw */
idx = loc(Type=2);
y = j(ncol(idx),1);
call randgen(y, "Erlang", 2, shape);
x[idx] = y; 

call series(t, mixDist); /* plot pdf */
call histogram(x) rebin={0.1, 0.2};

If you need to create many random samples, you can be even more efficient because the number of elements for each the exponential draw is Q = Binomial(N,p) and the number for the Erlang draw is N - Q. So no need to use the LOC function.  I just did it so that the code above is closer to the code on the blog post.

[SWEETSAS]

Thank you very much, Rick!

I must say that you are not from this earth! You are a genius to have figured that the distribution is a mixture by merely looking at the information I provided. It is indeed a mixture of exponential an gamma(2). In fact, you got the missing probability exactly right.  

But this is simplest form of the distribution that we work with. It is the flexibility of this distribution that makes it appealing. We work with several generalization of this distribution. The generalization allows for flexibility in modeling; the generalizations cannot easily be expressed  as a mixture. For instance, there is a generalization in which the CDF is raised to a power, so that when the power=1 then we have the parent distribution. In such cases, one has to result to CDF or Lambert function.  That is the reason we want a SAS program that implements Lambert function and CDF that we can easily modify when we develop generalizations of this distribution. 

  Again, you are great thinker. And thanks for sharing your knowledge.

J

[Rick_SAS]

You still haven't said which branch of the function you need.

There is a lot of information in the literature about how to efficiently compute the Lambert W function. The paper byChapeau-Blondeau and Monir (2002) seems to have a nice fast direct method for the lower branch, which has applications to generating data from the long-tailed generalized Gaussian distribution.  For the principal branch, an initial guess as described in Boyd (1998) looks promising, 

followed by one or two iterations of Halley's method.

Good luck to you.

[SWEETSAS]

Thanks very much much.

I need the negative branch of Lambert function. Here is quantile function expressed in terms of the negative branch of the Lambert function.

Quantile function (Q):

Q(p|alpha)=-1-(1/alpha)-(1/alpha)W_₁[(1+alpha)*(p-1)*exp(-(1+alpha))];

0<p<1; Where W_₁ where denotes the negative branch of the Lambert W function

P is essentially a unif(0,1)

Thanks again!

J

[Rick_SAS]

Using your current software, how long does does it take to compute 1E7 random variates by calling the quantile function (which calls the Lambert functioN) with alpha=1?

[SWEETSAS]

execution time is 3.65 seconds; CPU time is 0.15 seconds

[Rick_SAS · Accepted solution]

SAS/IML 14.3 (released with SAS 9.4M5) contains a built-in implementation of the Lambert W function.