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? 
 A: Naïvely, what it looks like you have is a GLM, depending on the variance of $\epsilon$. Furthermore, I believe 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[[1]]
    b_hat <- par[[2]] 
    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[[1]]
    b_hat <- par[[2]] 
    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
[1] 9.001184 2.999822

[1] 8.999761 3.000036
[1] 9.000211 2.999969
[1] 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

