# Non-linear regression for Lambert W-function in stata

I need to estimate parameters $\alpha, \beta$ and $k$ in equation

$y_i = \alpha + \beta \times W(k x_i)+\epsilon_i$

where ($x_i$,$y_i$) is the data and $W(\cdot)$ is Lambert W-function. This is the function with the property $W(xe^x)=x$.

$Var(\epsilon_i)=f(x_i)$ where $f(\cdot)$ is unknown. $x_i$ is within $[0,10000]$ range.

There are 300K observations.

What would be the easiest way to solve the problem?

• It looks like there is a closed form expression for $W$ (en.wikipedia.org/wiki/Lambert_W_function#Derivative), so you could probably just choose your favorite numerical optimization technique. Nonlinear least squares, gradient descent, etc. But I'm curious about what someone with numerical optimization experience has to say here – shadowtalker Sep 5 '16 at 5:08
• What are you assuming about $\epsilon$? (e.g. is variance roughly constant or does it change with x in some way?) – Glen_b -Reinstate Monica Sep 5 '16 at 6:18
• Is your biggest problem getting good starting values? Or is there an implementation problem? – Glen_b -Reinstate Monica Sep 5 '16 at 7:01
• Implementation. In particular, an efficient way to compute $W(kx_i)$ – P.Escondido Sep 5 '16 at 7:06
• What language are you using? R has lamW package ( a fast Rcpp implementation), python has scipy special functions, GSL library has it, boost library will get a lambertw implementation soon. – Georg M. Goerg Jan 4 '17 at 11:11

Naïvely, what it looks like you have is a GLM, depending on the variance of $\epsilon$. Furthermore, I beleive you have set up an indeterminate problem as regardless of what $k$ is, it can be offset by $\beta$. You may be better off ignoring $k$ completely, in which case take your $x_i$ values and use existing software (like my lamW package in R) to calculate $W(x_i)$ (call it $z_i$). Now run your GLM on $y_i = \alpha + \beta z_i$ and select the variance function (normal, gamma, etc.) which gives you your best fit.

Or, assuming $E(\epsilon) = 0$ you can try brute force it with something like this:

library(lamW)
library(nloptr)
Regress <- function(par, x, y){
a_hat <- par[]
b_hat <- par[]
y_hat <- a_hat + b_hat * lambertW0(x)
sum((y - y_hat) ^ 2)
}

Fit <- nloptr(x0 = c(2, 2), eval_f = Regress, x = x, y = y,
opts = list(algorithm = "NLOPT_LN_SBPLX", maxeval = 1e5,
ftol_abs = 1e-9, ftol_rel = 1e-9, tol_rel = 1e-7))
Fit$solution for (i in seq_len(3)){ Fit <- nloptr(x0 = Fit$solution, eval_f = Regress, x = x, y = y,
opts = list(algorithm = "NLOPT_LN_SBPLX", maxeval = 1e5,
ftol_abs = 1e-9, ftol_rel = 1e-9, tol_rel = 1e-7))
print(Fit$solution) }  For example, running the code below returns the simulated values. set.seed(12) library(lamW) library(nloptr) a <- 9 b <- 3 x <- runif(3e5, 0, 10000) e <- rnorm(3e5, 0, sqrt(x)) # Make Var(e) depend on x y <- a + b * lambertW0(x) + 0 * e Regress <- function(par, x, y){ a_hat <- par[] b_hat <- par[] y_hat <- a_hat + b_hat * lambertW0(x) sum((y - y_hat) ^ 2) } Fit <- nloptr(x0 = c(2, 2), eval_f = Regress, x = x, y = y, opts = list(algorithm = "NLOPT_LN_SBPLX", maxeval = 1e5, ftol_abs = 1e-9, ftol_rel = 1e-9, tol_rel = 1e-7)) Fit$solution
for (i in seq_len(3)){
Fit <- nloptr(x0 = Fit$solution, eval_f = Regress, x = x, y = y, opts = list(algorithm = "NLOPT_LN_SBPLX", maxeval = 1e5, ftol_abs = 1e-9, ftol_rel = 1e-9, tol_rel = 1e-7)) print(Fit$solution)
}

Fit


Returns

> Fit$solution  9.001184 2.999822  8.999761 3.000036  9.000211 2.999969  9.000082 2.999989 > Fit Call: nloptr(x0 = Fit$solution, eval_f = Regress, opts = list(algorithm = "NLOPT_LN_SBPLX",
maxeval = 1e+05, ftol_abs = 1e-09, ftol_rel = 1e-09, tol_rel = 1e-07),     x = x, y = y)

Minimization using NLopt version 2.4.2

NLopt solver status: 4 ( NLOPT_XTOL_REACHED: Optimization stopped because xtol_rel or xtol_abs (above) was reached. )

Number of Iterations....: 149
Termination conditions:  maxeval: 1e+05 ftol_abs: 1e-09 ftol_rel: 1e-09
Number of inequality constraints:  0
Number of equality constraints:    0
Optimal value of objective function:  5.39067587198454e-05
Optimal value of controls: 9.000082 2.999989