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?
 A: 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

A: 
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. 
