I have a question on standard errors for non-linear least squares in R. With the built-in function NLS and a hand-made function I get different SE and I don't understand why. I will try to expose below a simpler and replicable example of my original problem.
I want to estimate the mean $\mu$ and the standard deviation $\sigma$ of a normal cumulative distribution function.
Let's simulate the data for this example in R:
set.seed(12)
mu <- 5
sigma <- 2
x <- runif(100, min = 0, max = 10)
y <- pnorm(x + rnorm(1), mu, sigma)
df <- data.frame(y,var1)
I assume a true value for $\mu$ of 5 and 2 for $\sigma$. I add a normal noise to y which is a cdf of x.
I first fit a non-linear least squares with the built-in function NLS:
res1 <- nls(y ~ pnorm(x, mean = mu, sd = sigma),
data = df, start = c(mu = 1,sigma = 1), algorithm = "port")
summary(res1)
I get the following result that I was expected
Formula: y ~ pnorm(x, mean = mu, sd = sigma)
Parameters:
Estimate Std. Error t value
mu 5.04268491217399539295 0.00000000000000003347 150682292747528064
sigma 2.00000000000000000000 0.00000000000000004567 43795073347627272
Pr(>|t|)
mu <0.0000000000000002 ***
sigma <0.0000000000000002 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.00000000000000003995 on 98 degrees of freedom
Algorithm "port", convergence message: X-convergence (3)
Secondly, I tried my own LS estimation
f <- function(par){
mu = par[1]
sigma = par[2]
temp <- pnorm(df$x, mean = mu, sd = sigma)
sum1 <- sum((df$y-temp)^2)
return(sum1)
}
res2 <- optim(c(1,1),f,method = "L-BFGS-B",
lower=c(0,0), upper = c(100,100), hessian = TRUE)
print(res2)
I get similar result for the coefficients:
$par
[1] 5.042685 1.999999
$value
[1] 0.000000000001761186
$counts
function gradient
20 20
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
$hessian
[,1] [,2]
[1,] 2.8835336 -0.2273554
[2,] -0.2273554 1.5485147
But when I compute the Standard Errors, I don't get the same result with the NLS:
sqrt(diag(res2$hessian))
[1] 1.698097 1.244393
Does someone know why ?
Thank you !
Alexis