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.
5 ビュー (過去 30 日間)
古いコメントを表示
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.
3 件のコメント
Baptiste
2020 年 4 月 30 日
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
回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Get Started with Curve Fitting Toolbox についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!