How to correctly use lsqnonlin in conjunction with nlparci and nlpredci, by fitting to a weighted function and utilising the jacobian and residue output from lsqnonlin.

I have three (x,y) experimental data sets, representing data at three geographical locations, call them A [x_A, y_A] of size n_A, B [x_B, y_B] of size n_B and C [x_C, y_C] of size n_C, x_A~=x_B~=x_C, but all are contained within an interval [0,X_max], all y values are >0.

I have a model that describes this system, based on three parameters p,q and r. I have a function that takes a set of x values and returns y at each location A,B,C at the given x values. Call this function:

[ y_model_A, y_model_B, y_model_C ] = func( x, p, q, r )

I want to fit func by varying p,q and r to the three data sets. I have attempted this using lsqnonlin. To write the difference function to be used by lsqnonlin I did as followed

   function [ diff ] = diff_func( parameters )

          % grab parameters needed for func from 'parameters'
          p = parameters(1);
          q = parameters(2);
          r = parameters(3);

          % initialise diff vector
          diff = [];

          % solve analytically to find this solution over the timesteps found
          % in the experimental data, A,B and C.

          %%ret
          [ temp ] = func( x_A, p, q, r );
          y_model_A = temp(:,1);

          %%vit
          [ B_temp ] = func( x_B, p, q, r );
          y_model_B = temp(:,2);

          %%ret
          [ C_temp ] = func( x_C, p, q, r );
          y_model_C = temp(:,3);

          %%objective function evaluation
          diff = [ diff; objective_function( y_model_A, y_A ) ];
          diff = [ diff; objective_function( y_model_B, y_B ) ];    
          diff = [ diff; objective_function( y_model_C, y_C ) ];

end

Where I have defined my objective function as

function [ diff_partial ] = objective_function( data_model, data_exp )

data_model = log(data_model);
data_exp = log(data_exp);

% number
num = length(data_exp);

%%processing
if isempty(data_exp) == 0
    diff_partial = (1/num)*( data_model- data_exp )./data_exp;
end

end

I then evaluated lsqnonlin as such

    par_0 = [ p_0, q_0, r_0 ] % initial guesses for p,q and r 
    [ par, ~, residual, ~, ~, ~, jacobian ] = lsqnonlin( @diff_func, par_0 );

This works, I get nicely fitted p,q and r from the 'par' output. However I now want to use 'residual' and 'jacobian' to create confidence intervals for my parameters.

My question is, given what I have done here, how can I use nlparci and nlpredci in conjunction with 'residual' and 'jacobian' to produce confidence intervals for my parameters p,q and r (using nlparci) and for my function output (nlpredci).

I am concerned because both of these functions seem to prefer being used with nlinfit. I'm worried that because of how I have defined my difference function, specifically due to the weighting in my difference function that I will not be getting the intended result by utilizing these confidence interval functions with the generated jacobian and residues. I would use nlinfit but I need to constrain p,q and r to be >0, nlinfit gives me negative values.

Given the code I have written here how I can use nlparci and nlpredci and be sure that I am getting the correct (as in mathematically correct) answer? I would very much appreciate some help with this as I am very new to fitting data to models and have not used this MATLAB toolkit before.

EDIT: I have attached some code that demonstrates what I have described here for a function with artificial noise, if this function (despite it's definition of the objective function, with weights etc.) is giving me a statistically correct (approximately obviously) CI then I have correctly implemented these functions. Many thanks again in advance.

[jgg]

Have you considered bootstrapping your standard errors instead?

http://www.mathworks.com/help/stats/bootstrp.html

Basically, bootstrap like 500 solutions of your problem, and use the distribution of the solutions to calculate the CI for your problem.

You would avoid the problem of trying to get a statistically correct closed form, or trying to apply a CI method not designed for your problem. The downside is that it might be pretty slow; but it's a fairly robust technique which I think might work in your situation.

[Laurence hutton-smith]

This is very good advice, thank you very much for this link. I will look into this method!

[Baptiste]

I have exactly the same question than Laurence. I wheighted the error by my data because I had a decay with values on several decades and rigorous NLS fit could'nt fit properly

Error=(data_model-data_exp)./data.exp

I used lsqnonlin because I need bounds

Now I would like to use nlparci. Is it still working?

bootrstrap can be very slow since my function is slow and I have seven parameters.

What is your opinion?

Regards