Goodness of Fit - Error Bars

I'm currently using the fminsearch function in order to find a best fit curve to a non-analytic function. However, I would also like to be able to find error bars on my fit parameters, say at a 95% confidence level.

I don't have any experience with statistical analyses, and have been relying on built in fit functions to obtain error bars in the past. Is there any good way to do that without drastically altering my program?

To give a bit more information about my program: I am solving a system of equations with given parameters p(1), p(2), p(3), and then subtracting experimental data, x. The fminsearch function searches for optimal p by referencing a function (here, "functionname") which computes the quadratic mean (RMS) RMS(f(p)-x)

[poptimal, fval] = fminsearch(@functionname, p0);

I'd like to comment ahead of time that I don't want errorbars on the graph of my function or anything like that. I just want to know, for example, that p(1) is equal to 3e-6 +/- 2e-7 with 95% confidence.

Answer 1

[the cyclist]

I could be wrong, but I don't think that fminsearch() will be able to give you a canned goodness-of-fit error bar, because in a sense it does not "know" that you are doing a fit at all. You should theoretically be able to calculate it, but MATLAB isn't going to do it for you.

I believe the more standard approach for what you are doing would be to use the nlinfit() function in the Statistics Toolbox. That function produces error estimate output for the fitted parameters.

[Clement Wong]

I just looked into the nlinfit funciton, and I'm not sure I understand how it works. Currently, my function is as follows:

HanleCurve(B,Rch,L,D,tau,Rc,W,I) calculates a curve for all values of vector B, at given values for Rch,L,Rc,W, and I. The two parameters I am trying to fit are D and tau, as well as a amplitude multiplier, k, on HanleCurve (see explanation of HanleRMS below). However, this function is not analytic, and is arrived at by solving a 12x12 matrix for it's eigenvalues individually for each point in B.

HanleFit takes in a vector x0 such that x(1) = k0, x(2) = D0, x(3) = tau0 are the initial values from which to search. HanleFit has a subfunction, HanleRMS, which computes RMS(y-x(1)*HanleCurve(...,x(2),x(3)...)), where RMS computes the quadratic mean, and y is the experimental data. I have neglected to write out the known parameters of HanleCurve for simplicity.

HanleFit then calls the function fminsearch(@HanleRMS,x).

So, if I were to use nlinfit, I would replace fminsearch with nlinfit. However, I don't understand what the parameter X in nlinfit is supposed to represent, or how to write a function in the form: yhat = modelfun(b,X)

Can you explain what it means by X is a n-by-p matrix of p "predictors" for each point n? What is a predictor, and how do I find them for my system? Is this approach fine for this problem?

Answer 2

[Clement Wong]

After working on it for a while longer, I seem to have figured out how to use nlinfit. My data does indeed fit properly. Now my next question is whether or not I can get confidence intervals that take into account dependencies between my variables. A simple example is fitting y = (A+B)*x. I should get, for any data, that the confidence intervals for A and B are -inf to inf.

However, this function doesn't seem to do that. Any ideas?

[the cyclist]

I disagree with you. Here is a simple set of (x,y) data. In the first call to nlinfit(), I use the model y = P*x. In the second call, I used y = (P1+P2)*x. The first call works fine. The second call gives me a warning that the Jacobian is ill-conditioned, and the covariance matrix of the parameter estimates is of order 1e27. This is exactly the behavior I would expect. [Sorry that the code will probably be very poorly formatted, but there is no markup in comments.]

% Pretend data

x = 1:8;

y = 3*x + 0.4*rand(1,8);

% Well parameterized model

f = @(p,x) p(1)*x;

beta0 = [1];

[beta,r,j,covb,mse] = nlinfit(x,y,f,beta0)

% Poorly parameterized model

f = @(p,x) (p(1)+p(2))*x;

beta0 = [1 1];

[beta,r,j,covb,mse] = nlinfit(x,y,f,beta0)

[Clement Wong]

When I use the data set x = 1:5, y = 1:5, and the function f = @(p,x) (p(1)+p(2))*x, with beta0 = [1 1], I get a covariance matrix with values ranging only from -8.5 to 8.5. When I then run nlparci(beta,r,'covar',covb), I get confidence intervals only from -8.7 to 9.7. With such a poorly parameterized function, I would expect to get something astronomically huge, or, if MATLAB had a function to support it, infinity itself. Can you elaborate on what I'm doing wrong?