Re: confidence region and leasqr()

[Olaf Till]

On Wed, Apr 10, 2013 at 09:19:17PM +0400, Dmitry Shkirmanov wrote:
> Thanks for reply. I am not expert in this field. By google i found
> delta method
> ( 
> http://stats.stackexchange.com/questions/15423/how-to-compute-prediction-bands-for-non-linear-regression
> )
> and 
> <http://stats.stackexchange.com/questions/9833/constructing-95-confidence-interval-based-on-profile-likelihood>profile
> likelihood approach
> ( 
> http://stats.stackexchange.com/questions/9833/constructing-95-confidence-interval-based-on-profile-likelihood
> )
> 
> 
> Also, there is similar question in R mailing list archive (
> https://stat.ethz.ch/pipermail/r-help/2010-August/247918.html )

These references use the Delta method (which I did not know before) to
compute the variance of the computed model function values. I still
have not tried to verify or seen a proof that this is correct, but
probably it is.

It would actually not be correct to use this variance to compute
confidence intervals with the normal distribution, since it is not the
true variance, but itself a guess. So one has to use some other
statistic with a different distribution, but I don't know which.

But as long as the latter is not clarified, if you are willing to use
this way although it is not (quite) correct (your confidence intervals
will be too small), you could do the following:

- obtain 'covp', the covariance matrix of the parameters, and 'dfdp',
  the jacobian of the computed model function values against the
  parameters. You could use leasqr to get 'covp' and possibly the 'dfdp'
  function to get 'dfdp', but more straightforward is to use
  'nonlin_curvefit' and 'curvefit_stat' (note that these need the
  arguments of the user function in reversed order compared to leasqr;
  and the first line of the helptext of curvefit_stat contains the wrong
  function name as I've just seen, but otherwise should be correct).
Example:

octave:15> x = 0:5;
octave:16> y = [10, 5, 3, 2, 1.5, 1];
octave:17> f = @ (p, x) p(1) * exp (p(2) * x);
octave:18> pin = [1; 1];
octave:19> [p, fy, cvg, outp] = nonlin_curvefit (f, pin, x, y);
octave:20> info = curvefit_stat (f, p, x, y, optimset ("ret_dfdp", true, 
"ret_covp", true, "objf", "wls"));

- then get the variances as:

octave:24> v = diag (info.dfdp * info.covp * info.dfdp.');

- choose your level of significance:

octave:25> p = .05;

- get half width of confidence intervals: (as stated above, not quite
  correct):

octave:26> d = norminv (1 - p / 2, 0, sqrt (v));

- compute upper and lower borders of your confidence bands:

octave:30> upper = fy(:) + d;
octave:31> lower = fy(:) - d;

- and plot with something like:

octave:32> plot (x, lower, x, fy, x, upper)

I'd still like to know which distribution should be used instead, then
I probably could automatize this in residmin_stat and curvefit_stat...

Olaf

-- 
public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net

signature.asc
Description: Digital signature