2
$\begingroup$

I am using implementations of the Levenberg-Marquardt algorithm for non-linear least squares regression based on MINPACK-1 utilizing either the R function nlsLM() from minpack.lm or an implementation in C using the mpfit library.

I have observed that while I come up with comparable parameter estimates in both cases, the calculated uncertainty is drastically different in the two approaches. Below you find the code of the two approaches:

In R:

library(minpack.lm)

x <- seq(1,9)
y <- c(2.93,1.79,1.03,0.749,0.562,0.218,0.068,0.155,0.03)
table <- data.frame(x,y)

fit <- nlsLM(y ~ N * exp( -rate * x ), data = table, start = list(rate = 0.1, N = 1)

summary(fit)

which returns the estimates and standard errors:

rate    0.47825, std.Err: 0.02117
   N    4.69723, std.Err: 0.18950

In C:

#include <math.h>
#include <stdio.h>
#include <string.h>

#include "cmpfit-1.2/mpfit.h"

// Data structure to keep my vectors
typedef struct my_data {

  double *x;
  double *y;
  int n;

} my_data;

// defining the function that is to be fitted
int my_exp(int m, int n, double *p, double *dev, double **deriv, my_data *data)
{

  int i;
  double *x, *y, f;

  x = data->x;
  y = data->y;

  for (i = 0; i < m; i++) {

    f = p[0] * exp(-p[1]*x[i]);
    dev[i] = y[i] - f;

  }

  return 0;
}


int main(void)
{

  // define test data
  my_data data;
  data.x = (double[]){1,2,3,4,5,6,7,8,9};
  data.y = (double[]){2.93,1.79,1.03,0.749,0.562,0.218,0.068,0.155,0.032};
  data.n = 9;

  // start values
  double p[2] = {1.0, 0.1};

  mp_result result;
  memset(&result, 0, sizeof(result));

  double perr[2];
  result.xerror = perr;

  int ret;

  ret = mpfit((mp_func)my_exp, data.n, 2, p, 0,0, &data, &result);


  printf("rate:   %f, Error: %f\n", p[1], perr[1]); 
  printf("   N:   %f, Error: %f\n", p[0], perr[0]);

  return 0;
}

which returns:

rate:   0.479721, Error: 0.270632
   N:   4.707439, Error: 2.422189

As you can see, the errors returned by mpfit are about 12x higher then from R. I also implemented the same thing using GSL which produces results identical to mpfit. Unfortunately I am not very comfortable with the R code base, so it would be great if anyone has any inside into what nlsLM in R does differently then mpfit and GSL when calculating the errors.

* EDIT 1 *

So, I've made a few discoveries I'd like to share and hopefully get some input on. The errors returned by mpfit in the C implementation actually correspond to what one would get in R using:

diag(sqrt(summary(fit)$cov.unscaled))

i.e. the square root of the diagonal of the unscaled covariance matrix. The standard error estimates for the parameters actually returned by summary(fit) are apparently scaled by the residual standard error of the regression, so that:

sqrt(diag(summary(fit)$cov.unscaled) * s$sigma^2)

actually returns the values also observed in summary(fit).

Why is this scaling done? What is the statistical reason behind it? Any input would be very welcome!

* EDIT 2 *

So, I think I figured it out now and I'll put it here for future reference. What is returned by the C implementation as error, as well as stored in cov.unscaled in R is just the equivalent of solve(t(J) %*% J) which needs to be scaled by sigma^2 to get the variance of the parameter estimates. Thanks @dave for the help!

$\endgroup$
1
$\begingroup$

I fit your data with AD Model Builder. The R code has the right std devs (almost) but poor parameter estimates. The C code has good parameter estimates but the wrong std devs You need to multiply by sqrt(9/7) to convert from fisher information to nls stdevs

index name   value      std.dev 
  1     N   4.7079e+00 1.7120e-01 
  2    rate 4.7979e-01 1.9154e-02  
$\endgroup$
  • $\begingroup$ Hi @dave, thank you very much for your response. Could you clarify what I need to multiply by sqrt(9/7)? Simply multiplying the error from the C code clearly doesn't yield the errors from R or your example, as sqrt(9/7) > 1. Also, do you have any inside on how the standard error in R (or AD model builder for that matter) is actually calculated compared to the one mpfit puts out? Thanks! $\endgroup$ – clw Dec 5 '16 at 16:42
  • $\begingroup$ please see my edits in the original question. Any comment on that would be very welcome. $\endgroup$ – clw Dec 5 '16 at 17:20
  • $\begingroup$ The AD Model Builder estimates for the std devs are based on the inverse of the Hessian of the log-likelihood function. For nls the log-likelihood has special structure involving the MLE for the variance of the errors. This is the sum of the squared residuals divided by the number of observations in this case 9. For NLS this is corrected by replacing the number of observations by the number of observations 9 minus the number of parameters 2 =7. So to get that you multiply by 9/7 and for std dev take square root. So you multiply the AD Model Builder std devs by sqrt(9/7) to get the nls ones. $\endgroup$ – dave fournier Dec 5 '16 at 17:34
  • $\begingroup$ Ah, I see. Thats to compare between AD Model Builder and the nls in R. I was trying to figure out how the errors reported by R compare to those reported by the mpfit in C. I have made some discoveries on that front which I included in the original question as an edit. It would be great if you could comment on that too. Thanks! $\endgroup$ – clw Dec 5 '16 at 17:41

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.