Re: confidence region and leasqr()

[Olaf Till]

On Fri, Apr 12, 2013 at 10:43:10PM +0400, Dmitry Shkirmanov wrote:
> 
> >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.
> If i understood correctly the Delta method is a bad approach.

It's not a bad approach; it seems to provide a fist approximation of
the (co-)variance of computed model function values, using the
covariance matrix of the computed parameters. But of the latter you
usually have only a guess.

After reading a bit on this, I came up with this:

- example data and model function as previously:

octave:1> x = 0:5;
octave:2> y = [10, 5, 3, 2, 1.5, 1];
octave:3> f = @ (p, x) p(1) * exp (p(2) * x);
octave:4> pin = [1; 1];

- get a fit:

octave:5> [p, fy, cvg, outp] = nonlin_curvefit (f, pin, x, y);

- obtain covariance matrices of parameters ('covp') and of data
  ('covd'), and jacobian of model function against parameters
  ('dfdp'):

octave:6> info = curvefit_stat (f, p, x, y, optimset ("ret_dfdp", true, 
"ret_covp", true, "ret_covd", true, "objf", "wls"));
octave:8> dfdp = info.dfdp; covp = info.covp; covd = info.covd;

- get guessed variances of all points of the computed model function:

octave:9> s2 = diag (dfdp * covp * dfdp.');

- since you want prediction intervals, not confidence intervals, add
  the guessed variance of data:

octave:10> s2 += diag (covd);

- Now (residuals.^2 ./ s2) should be distributed as
  Chi-square(1). Since the used procedures assumed covd up to a factor
  which is guessed, and computed covp from covd, get rid of the
  guessed factor by computing (assuming near linearity of the
  parameters):

octave:11> resid = (fy - y)(:);
octave:12> m = 6; # number of data points
octave:13> n = 2; # number of parameters
octave:15> Z = (m - n) ./ (s2 * (resid.' * covd * resid));

- residuals.^2 * Z should now be distributed as F(1, m - n). So
  compute half width of confidence band as:

octave:30> hw = sqrt (finv (1 - alpha, 1, m - n) ./ Z);

- and again, plot with something like:

octave:32> plot (x, fy(:) - hw, x, fy(:), x, fy(:) + hw);

Note that this is a _pointwise_ confidence band, so while you can make
a statement for each point, you _cannot_ say all observations of a
future experiment will be within this band with this level of
confidence. For the latter, you'd need a simultaneous confidence
band. While there exist methods to obtain these, I have not studied
them (and probably will not); in any case it seems to me that,
counterintuitively, the width of simultaneous confidence bands will
increase with the number of sample points.

Note also that there are assumptions in computing the bands which are
not necessarily fullfilled (near linearity (see above), covariance
matrix of data defined up to a factor by weights given by user).

Olaf

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

signature.asc
Description: Digital signature