Lambert function

[Nico Gatti]

Hi everyone,
I would need to compute the Lambert function W(a) in Stata, but I cannot find a command to do this.
Is anybody aware of a fast way to compute this value?
Thanks a lot,
N.

[Leonardo Guizzetti]

This is beyond my capability, but Mata provides facility for computations using complex numbers.

Example:

mata
  eltype(5 + 3i)
  (1i)^2
end

[Nico Gatti]

Thanks a lot!
Actually, I have never used Mata.
Estimating the Lambert function means solving the equation x*e^x = a, where a is known.
Unfortunately, I am not able to do this in Stata.

[Hua Peng (StataCorp)]

If you are using Stata 16 and Python is installed on the machine, the simplest way to compute lambert W is:

python:
from scipy.special import lambertw

w = lambertw(1)
w

w = lambertw(1+5j)
w
end

Running it in Stata:

. python:
----------------------------------------------- python (type end to exit) -----
>>> from scipy.special import lambertw
>>>
>>> w = lambertw(1)
>>> w
(0.5671432904097838+0j)
>>>
>>> w = lambertw(1+5j)
>>> w
(1.2399002298584543+0.8002535801631443j)
>>> end
-------------------------------------------------------------------------------

[Nico Gatti]

Hua Peng (StataCorp) Thank you so much!!!
It would be great, but at the moment my institution has available only Stata 15... Isn't there a solution which could be used in this case?Thank you so much

[Hua Peng (StataCorp)]

The following is a Mata port of TOMS433 by John Burkardt, https://people.sc.fsu.edu/~jburkardt...s443/toms443.c

mata:


/* 
    The code is a straight Mata port of  the C version of TOMS443
    
    Original Author:

        John Burkardt

    Reference:

    Brian Hayes,
        "Why W?",
        The American Scientist,
        Volume 93, March-April 2005, pages 104-108.

    Eric Weisstein,
        "Lambert's W-Function",
        CRC Concise Encyclopedia of Mathematics,
        CRC Press, 1998.

    Stephen Wolfram,
        The Mathematica Book,
        Fourth Edition,
        Cambridge University Press, 1999,
        ISBN: 0-521-64314-7,
        LC: QA76.95.W65.
    
    John Burkardt,
        https://people.sc.fsu.edu/~jburkardt/c_src/toms443/toms443.html
        https://people.sc.fsu.edu/~jburkardt/c_src/toms443/toms443.c
    
    Licensing:    
        The original code is distributed under the GNU LGPL license.
        Hence the Mata code is also distributed under the GNU LGPL license.
*/

real scalar wew_a (real scalar x, real scalar en )
{
  real scalar f
  real scalar temp
  real scalar temp2
  real scalar wn
  real scalar y
  real scalar zn
  real scalar c1, c2, c3, c4

  c1 = 4.0 / 3.0;
  c2 = 7.0 / 3.0;
  c3 = 5.0 / 6.0;
  c4 = 2.0 / 3.0;

    // Initial guess
  f = log ( x );

  if ( x <= 6.46 ) {
    wn = x * ( 1.0 + c1 * x ) / ( 1.0 + x * ( c2 + c3 * x ) )
    zn = f - wn - log (wn)
  }
  else {
    wn = f
    zn = - log(wn)
  }

    // Iteration 1
  temp = 1.0 + wn
  y = 2.0*temp*(temp + c4*zn) - zn
  wn = wn*(1.0 + zn*y/(temp*(y-zn)))
  
    // Iteration 2.
  zn = f - wn - log(wn)
  temp = 1.0 + wn
  temp2 = temp + c4*zn
  
  en = zn * temp2 / (temp * temp2 - 0.5 * zn )
  wn = wn * ( 1.0 + en )
  return(wn)
}

end

The function only computes the principle branch for real values. The following are some results:

mata:

: w = wew_a(1, en=.)

: w
  .5671432904

: w = wew_a(2, en=.)

: w
  .852605502

: w = wew_a(0.5, en=.)

: w
  .3517337112

: w = wew_a(0.75, en=.)

: w
  .4691502107

end