# Estimate $\lambda$ in the fitting line $x^\lambda$, where $x \in [0, 1]$

## Problem

I would like to estimate $$\lambda$$ in the fitted line $$x^\lambda$$, where $$x \in [0, 1]$$.

Note that the following R code generates "concave" growth as x increases from $$0$$ to $$1$$.

Lambda = 1/2.42

x = rbeta(1e4, shape1=2,shape2=2)
y = x^Lambda + rnorm(1e4, sd=.1)

plot(x,y)


## Try

My approach is to build a loss function, the L2 loss to find the optimal $$\lambda$$, like the following R code.

SumSq = function(lam) sum((y - x^lam )^2)
optimize(SumSq, c(0,1), tol=1e-4, maximum = F )


This code assumes I know $$\lambda \in (0,1)$$, and optimize gives me a quite accurate result.

But this approach cannot give any statistical asymptotic results such as CI, which is useful information if I would like to consider the uncertainty.

Is there any standard way to do it?

• i really dont know why someone would down vote and not leave a comment. This is a reasonable question with a minimally reproducible example. Commented Mar 7, 2019 at 0:48
• @forecaster I agree and found your answer to this question to be edifying, so I am similarly puzzled. Commented Mar 7, 2019 at 2:52
• The natural approach to estimations with this problem would involve calling a nonlinear least squares routine (nls in R - there are a number of somewhat more robust alternatves, but this is usually sufficient). Commented Mar 7, 2019 at 5:25
• When the model is expressed in the form $Y= \exp(\lambda \log x) + \varepsilon,$ this question is seen to be the same as many other questions asked (and answered) on this site about nonlinear least-squares fitting: this gives you an additional set of resources to consult.
– whuber
Commented Mar 7, 2019 at 16:29
• "But this approach cannot give any statistical asymptotic results such as CI, which is useful information if I would like to consider the uncertainty." Is this the core of the question? This is not so clear from the title which seems to point to point estimates (pun intended). Commented Mar 7, 2019 at 19:43

Firstly I would suggest you read an excellent book by Ben Bolker entitled Ecological Models and Data in R. Written for ecologists , I think this is one of the best books on practical data analysis regardless of background, I'm an engineer and I have used it a lot. I'm sure there are other mathematical statistics book, however this book is by far the most practical book that I have read. He also has a package BBMLE which you might want to check. The book unlike any other goes in to maximum likelihood estimation, profile likelihood and confidence interval estimation and all that.

Coming back to your problem, you need to write a log-likelihood of the function that you are trying to fit the data. Obviously you are assuming the error to be normally distributed as in the simulation, so the Log likelihood is:

$$Log\ Likelihood (Lambda,sd) = -\frac{n}{2} log(sd) - \frac{n}{2} log(2\pi)- \frac{(y-x^{Lambda})^2}{2sd}$$

You need to maximize the above equation using an optimization routine such as optim in R. Use the Hessian matrix from the optimization to assess the uncertainty of your estimated log likelihood function i.e., how steep or how flat the curvature of your function is at the optimal point. If the function is steep which implies less uncertainties at optimal point you would have a tighter confidence band on your parameter estimates, on the other hand if its flat you would have a wider confidence interval. Hessian, Fisher Information matrix would help you calculate the standard error and confidence interval.

Here is how you do it in R:

    set.seed(8345)

Lambda = 1/2.42

x = rbeta(1e4, shape1=2,shape2=2)
y = x^Lambda + rnorm(1e4,mean = 0, sd=.1)

plot(x,y)

## Write Log Likelihood function

log.lik <- function(theta,y,x){

Lam <- theta[1]
sigma2 <- theta[2]

# sample size
n <-  length(y)

#error
e<-y-(x^Lam)

#log likelihood
logl<- -.5*n*log(2*pi)-.5*n*log(sigma2)-((t(e)%*%e)/(2*sigma2))

return(-logl) # R optim does minimize so to maximize  multiply by -1
}

## Estimate Paramters thru maximum likelihood

max.lik <- optim(c(1,1), fn=log.lik, method = "L-BFGS-B", lower = c(0.00001,0.00001), hessian = T,y=y,x=x)

# Lambda
Lam <- max.lik$par[1] #0.4107119 #Fisher Information MAtrix fisher_info<-solve(max.lik$hessian)
prop_sigma<-sqrt(diag(fisher_info))

## Estimate 95% Confidence Interval
upper<-max.lik$$par+1.96*prop_sigma lower<-max.lik$$par-1.96*prop_sigma

interval<-data.frame(Parameter = c("Lambda","sd"),value=max.lik$par, lower=lower, upper=upper) interval  • This is a generalized linear model so it seems more obvious to use the glm function and avoid a general non-linear optimizer with numerical differentiation. Commented Mar 7, 2019 at 21:10 Is there any standard way to do it? If you think that following is a good approximation (the model you simulate from as an example) $$y_i = x_i^\lambda +\epsilon_i, \qquad \epsilon_i\sim N(0,\sigma^2)$$ i.e., $$\log(E(y_i))=\lambda \log x_i$$ then you can use glm with family = gaussian("log") which is exactly this model lambda <- 1/2.42 set.seed(49564503) n <- 1e4 x <- rbeta(n, shape1 = 2,shape2 = 2) y <- x^lambda + rnorm(n, sd = .1) fit <- glm(y ~ log(x) - 1, family = gaussian("log"), start = c(0, 1)) summary(fit) #R #R Call: #R glm(formula = y ~ log(x) - 1, family = gaussian("log"), start = 1) #R #R Deviance Residuals: #R Min 1Q Median 3Q Max #R -0.42452 -0.06767 0.00090 0.07020 0.34418 #R #R Coefficients: #R Estimate Std. Error t value Pr(>|t|) #R log(x) 0.410666 0.001807 227.3 <2e-16 *** #R --- #R Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #R #R (Dispersion parameter for gaussian family taken to be 0.01033428) #R #R Null deviance: 1057.81 on 10000 degrees of freedom #R Residual deviance: 103.33 on 9999 degrees of freedom #R AIC: -17341 #R #R Number of Fisher Scoring iterations: 5  Notice that both the dispersion parameter and coefficient estimate match as expected. Now, you can make confidence interval with the standard error above or by using confint. • +1 There seem to be a - 1 missing in the model formula although it's there in the output from summary(fit). Commented Mar 7, 2019 at 21:51 • Thanks. I forgot to add that. I would prefer to keep the intercept but then the model would be more general than the model the OP simulates from. Commented Mar 7, 2019 at 22:17 • Your analysis confuses the model$E[\log y] = \lambda \log x$with the different model$\log(E[y]) = \lambda \log x.$When the variance of$\epsilon$is much smaller than$x^{2\lambda}$it's probably an OK approximation, but otherwise it's important to pay attention to the distinction--and to offer a different solution. – whuber Commented Mar 8, 2019 at 0:28 • @whuber doesn't the log link function express the following? $$\log(E[y]) = X\beta$$ Then using$X = x^\prime = \log(x)$you get: $$\log(E[y]) = \beta x^\prime \quad \rightarrow \quad E[y] = e^{ \beta x^\prime} = e^{ \beta \log(x)} = x^\beta$$ Commented Mar 8, 2019 at 0:41 • @whuber, I do not see why$E [\log (y)] $matters. The OP is using the model $$y|x \sim N (x^\lambda, \sigma^2)$$ where$\mu = x^\lambda = \exp ( \lambda \log (x))\$. So $$\log (E [y|x]) = \log (\mu) = \lambda \log (x)$$ Commented Mar 8, 2019 at 7:23