3
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – shadowtalker Sep 5 '16 at 5:08
  • $\begingroup$ What are you assuming about $\epsilon$? (e.g. is variance roughly constant or does it change with x in some way?) $\endgroup$ – Glen_b -Reinstate Monica Sep 5 '16 at 6:18
  • $\begingroup$ Is your biggest problem getting good starting values? Or is there an implementation problem? $\endgroup$ – Glen_b -Reinstate Monica Sep 5 '16 at 7:01
  • 2
    $\begingroup$ Implementation. In particular, an efficient way to compute $W(kx_i)$ $\endgroup$ – P.Escondido Sep 5 '16 at 7:06
  • 1
    $\begingroup$ 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. $\endgroup$ – Georg M. Goerg Jan 4 '17 at 11:11
1
$\begingroup$

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[[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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.