Unable to data fitting with Lsqnonlin and estimate the coefficient. Need help fitting the data.

Question

Anand Ra 2021-6-26

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/865955-unable-to-data-fitting-with-lsqnonlin-and-estimate-the-coefficient-need-help-fitting-the-data

评论： Matt J 2021-7-2

I am unable to get a good fit. I have to keep guessing the initial assumption d0 to manually get the better fit. Can I please get some help on what is wrog with my code. Btw, I am trying to determine the best coefficient by fitting the y_obs data array into an arbitrary function. Should I add constrains? If so, how?Thank you.

t1 = [0:300:21600]';
y_obs = [         
    0
   1.6216e-01
   2.9813e-01
   4.0805e-01
   4.9338e-01
   5.5928e-01
   6.0894e-01
   6.5506e-01
   6.8876e-01
   7.1773e-01
   7.4284e-01
   7.6433e-01
   7.8448e-01
   7.9912e-01
   8.1484e-01
   8.2950e-01
   8.3760e-01
   8.4707e-01
   8.5767e-01
   8.6299e-01
   8.6862e-01
   8.7433e-01
   8.7879e-01
   8.8736e-01
   8.9152e-01
   8.9903e-01
   9.0343e-01
   9.0984e-01
   9.1344e-01
   9.1615e-01
   9.2046e-01
   9.2313e-01
   9.2640e-01
   9.2992e-01
   9.3134e-01
   9.3242e-01
   9.3616e-01
   9.3986e-01
   9.4201e-01
   9.4145e-01
   9.4434e-01
   9.4252e-01
   9.4131e-01
   9.4249e-01
   9.4283e-01
   9.4395e-01
   9.4355e-01
   9.4690e-01
   9.4919e-01
   9.5255e-01
   9.5626e-01
   9.6282e-01
   9.6283e-01
   9.6672e-01
   9.6439e-01
   9.6887e-01
   9.7099e-01
   9.7529e-01
   9.8034e-01
   9.8302e-01
   9.8494e-01
   9.8828e-01
   9.8769e-01
   9.9130e-01
   9.9108e-01
   9.9422e-01
   9.9561e-01
   9.9737e-01
   1.0002e+00
   1.0034e+00
   1.0044e+00
   1.00329
   1.0];
d0 = 1.8055e-10;
fun = @(d) ypred(d, t1) - y_obs;
best_d = lsqnonlin(fun, d0)
predicted_y = ypred(best_d, t1);
plot(t1, y_obs, '*', t1, predicted_y, '-');
legend({'observed', 'predicted'})
title('lsqnonlin data fitting')
function y_pred = ypred(d, t1)
    a=0.0011;
    gama = 0.01005;
    
    L2 = zeros(14,1);
    L3 = zeros(100,1);
    L4 = zeros(100,1);
    L5 = zeros(100,1);
    S= zeros(73,1);
    y_pred = zeros(73,1);
    
    % t = 0; 
    
    L1 = ((8*gama)/((pi*(1-exp(-2*gama*a)))));
    format longE
    
    k =1;
    for t = t1(:).'
        for n=0:1:100
            L2(n+1) = exp((((2*n + 1)^2)*-d*pi*pi*t)/(4*a*a));
            L3(n+1) = (((-1)^n)*2*gama)+(((2*n+1)*pi)*exp(-2*gama*a))/(2*a);
            L4(n+1)= ((2*n)+1)*((4*gama*gama)+((((2*n)+1)*pi)/(2*a))^2);
            L5(n+1) = ((L2(n+1)*L3(n+1))/L4(n+1));
        end
        S((t/300) +1) = sum(L5);
        y_pred((t/300)+1)= 1 -(L1*S((t/300) +1)); % predicted data 
    end
end

2 个评论
显示无隐藏无

Matt J 2021-6-26

在 MATLAB Online 中打开

Data.mat

Looks like a pretty good fit to me. Also, it is easy to verify via plots that the minimum was coorectly found

load Data

d0 = 1.8055e-10;

fun = @(d) ypred(d, t1) - y_obs;

best_d = lsqnonlin(fun, d0)

Local minimum possible. lsqnonlin stopped because the size of the current step is less than the value of the step size tolerance.

best_d =

1.566238696858825e-10

pfun=@(d) norm(fun(d));

fplot(pfun,[0,5e-10])

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

hold on; plot(best_d*[1,1],ylim,'--rx'); hold off

Anand Ra 2021-6-26

Thanks. However, for instance if I give the initial assumption d0 = 3.2e-10, I get a perfect fit like the below plot. The issue is that I have different data sets (y_obs) and its not easy to keep guessing the intial assumption every single time to get that visually good fit. Which is why asking if I can do something like add constrains?

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2021-6-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/865955-unable-to-data-fitting-with-lsqnonlin-and-estimate-the-coefficient-need-help-fitting-the-data#answer_734275

Since it's a fairly inexpensive 1D function , I think you should just be able to do a 1D grid search over a suitable interval to find a reasonable d0.

16 个评论
显示 14更早的评论隐藏 14更早的评论

Anand Ra 2021-6-27

在 MATLAB Online 中打开

Thank you Matt. That definitely helps with the initial guess. However, wondering why I am getting even better fit if I give d0 = 3.2e-10 with lsqnonlin model. Btw, below is the update code:

t1 = [0:300:21600]';
y_obs = [         
    0
   1.6216e-01
   2.9813e-01
   4.0805e-01
   4.9338e-01
   5.5928e-01
   6.0894e-01
   6.5506e-01
   6.8876e-01
   7.1773e-01
   7.4284e-01
   7.6433e-01
   7.8448e-01
   7.9912e-01
   8.1484e-01
   8.2950e-01
   8.3760e-01
   8.4707e-01
   8.5767e-01
   8.6299e-01
   8.6862e-01
   8.7433e-01
   8.7879e-01
   8.8736e-01
   8.9152e-01
   8.9903e-01
   9.0343e-01
   9.0984e-01
   9.1344e-01
   9.1615e-01
   9.2046e-01
   9.2313e-01
   9.2640e-01
   9.2992e-01
   9.3134e-01
   9.3242e-01
   9.3616e-01
   9.3986e-01
   9.4201e-01
   9.4145e-01
   9.4434e-01
   9.4252e-01
   9.4131e-01
   9.4249e-01
   9.4283e-01
   9.4395e-01
   9.4355e-01
   9.4690e-01
   9.4919e-01
   9.5255e-01
   9.5626e-01
   9.6282e-01
   9.6283e-01
   9.6672e-01
   9.6439e-01
   9.6887e-01
   9.7099e-01
   9.7529e-01
   9.8034e-01
   9.8302e-01
   9.8494e-01
   9.8828e-01
   9.8769e-01
   9.9130e-01
   9.9108e-01
   9.9422e-01
   9.9561e-01
   9.9737e-01
   1.0002e+00
   1.0034e+00
   1.0044e+00
   1.00329
   1.0];
%************
% Guessing the initial assumption d0 by finding the minimum using min function
%Implementing fminbnd instead of lsqnonlin
pfun = @(d) norm( ypred(d, t1) - y_obs);
dsamps=linspace(0,5e-10,50);
[~,imin]=min(  arrayfun(pfun,dsamps) );
[best_d,fval,exitflag]=fminbnd(pfun,dsamps(imin-1), dsamps(imin+1),...
                              optimset('TolX',1e-14));
exitflag,
best_d,
%****************
%Voiding lsqnonlin
%d0 = 1.8055e-10;
%fun = @(d) ypred(d, t1) - y_obs;
%best_d = lsqnonlin(fun, d0)
%************
predicted_y = ypred(best_d, t1);
plot(t1, y_obs, '*', t1, predicted_y, '-');
legend({'observed', 'predicted'})
title('lsqnonlin data fitting')
function y_pred = ypred(d, t1)
    a=0.0011;
    gama = 0.01005;
    
    L2 = zeros(14,1);
    L3 = zeros(100,1);
    L4 = zeros(100,1);
    L5 = zeros(100,1);
    S= zeros(73,1);
    y_pred = zeros(73,1);
    
    % t = 0; 
    
    L1 = ((8*gama)/((pi*(1-exp(-2*gama*a)))));
    format longE
    
  
    for t = t1(:).'
        for n=0:1:100
            L2(n+1) = exp((((2*n + 1)^2)*-d*pi*pi*t)/(4*a*a));
            L3(n+1) = (((-1)^n)*2*gama)+(((2*n+1)*pi)*exp(-2*gama*a))/(2*a);
            L4(n+1)= ((2*n)+1)*((4*gama*gama)+((((2*n)+1)*pi)/(2*a))^2);
            L5(n+1) = ((L2(n+1)*L3(n+1))/L4(n+1));
        end
        S((t/300) +1) = sum(L5);
        y_pred((t/300)+1)= 1 -(L1*S((t/300) +1)); % predicted data 
    end
end

Anand Ra 2021-6-27

在 MATLAB Online 中打开

My bad; I am very sorry, I think I am incorrect on that instance, as I couldnt reporduce the fit as I mentioned earlier. Not sure what was going on.

Btw, I tried with a diffferent data set ( y_obs). And below the best fit that I can obtain? Its barely fitting the points. Can it be constrained any further?

t1 = [0:300:25200]';
y_obs = [         
   0
0.0127027
0.0378204
0.0823202
0.135036
0.185492
0.235428
0.280011
0.324162
0.36309
0.398897
0.434106
0.464692
0.492385
0.519071
0.543569
0.568757
0.58709
0.606308
0.624904
0.644106
0.659982
0.673968
0.68761
0.701253
0.714016
0.72741
0.740307
0.752054
0.761136
0.77098
0.782929
0.794314
0.803362
0.812927
0.823867
0.832376
0.839823
0.847113
0.85366
0.859776
0.866798
0.872073
0.877704
0.882976
0.888673
0.893651
0.898524
0.903607
0.908458
0.912884
0.916185
0.919978
0.924837
0.929271
0.931731
0.934375
0.937633
0.941785
0.94324
0.944147
0.948995
0.952695
0.953771
0.956396
0.957976
0.961972
0.965774
0.968082
0.970771
0.97097
0.974649
0.977063
0.97848
0.979965
0.980716
0.983084
0.98617
0.988391
0.989767
0.991297
0.993934
0.996887
0.996628
1.0
];
%************
% Guessing the initial assumption d0 by finding the minimum using min function
%Implementing fminbnd instead of lsqnonlin
pfun = @(d) norm( ypred(d, t1) - y_obs);
dsamps=linspace(0,5e-10,50);
[~,imin]=min(  arrayfun(pfun,dsamps) );
[best_d,fval,exitflag]=fminbnd(pfun,dsamps(imin-1), dsamps(imin+1),...
                              optimset('TolX',1e-14));
exitflag,
best_d,
%****************
%Voiding lsqnonlin
%d0 = 1.8055e-10;
%fun = @(d) ypred(d, t1) - y_obs;
%best_d = lsqnonlin(fun, d0)
%************
predicted_y = ypred(best_d, t1);
plot(t1, y_obs, '*', t1, predicted_y, '-');
legend({'observed', 'predicted'})
title('fminbnd data fitting')
%title('lsqnonlin data fitting')
function y_pred = ypred(d, t1)
    a=0.0011;
    gama = 0.01005;
    
    L2 = zeros(14,1);
    L3 = zeros(100,1);
    L4 = zeros(100,1);
    L5 = zeros(100,1);
    S= zeros(85,1);
    y_pred = zeros(85,1);
    
    % t = 0; 
    
    L1 = ((8*gama)/((pi*(1-exp(-2*gama*a)))));
    format longE
    
  
    for t = t1(:).'
        for n=0:1:100
            L2(n+1) = exp((((2*n + 1)^2)*-d*pi*pi*t)/(4*a*a));
            L3(n+1) = (((-1)^n)*2*gama)+(((2*n+1)*pi)*exp(-2*gama*a))/(2*a);
            L4(n+1)= ((2*n)+1)*((4*gama*gama)+((((2*n)+1)*pi)/(2*a))^2);
            L5(n+1) = ((L2(n+1)*L3(n+1))/L4(n+1));
        end
        S((t/300) +1) = sum(L5);
        y_pred((t/300)+1)= 1 -(L1*S((t/300) +1)); % predicted data 
    end
end

Anand Ra 2021-7-2

在 MATLAB Online 中打开

Hello Matt, a follow up question. Another model seems to produce ( the author cannot share the details of the code/algorithm) a different fit for the same data set as mine. The fit I obtained using fminbnd and lsqnonlin per your guidance is 1) and 2). The fit obtained by the model used by author X is 3)

Looks like the fitting has been constrained by author in 3). Any thoughts on that?

1)

estimated d =1.57e-10 m^2/s

2)estimated d =1.62e-10 m^2/s

3) Author's fit

estimated d =1.9e-10 m^2/s

I observe the the predicted data fit is not reaching one, but its passing through the maximum number of data points.

% Non-linear data fitting to estimate the diffusion coefficent
% Solver = lsqnonlin
t1 = [0:300:21600]'; 
y_obs = [         
    0
   1.6216e-01
   2.9813e-01
   4.0805e-01
   4.9338e-01
   5.5928e-01
   6.0894e-01
   6.5506e-01
   6.8876e-01
   7.1773e-01
   7.4284e-01
   7.6433e-01
   7.8448e-01
   7.9912e-01
   8.1484e-01
   8.2950e-01
   8.3760e-01
   8.4707e-01
   8.5767e-01
   8.6299e-01
   8.6862e-01
   8.7433e-01
   8.7879e-01
   8.8736e-01
   8.9152e-01
   8.9903e-01
   9.0343e-01
   9.0984e-01
   9.1344e-01
   9.1615e-01
   9.2046e-01
   9.2313e-01
   9.2640e-01
   9.2992e-01
   9.3134e-01
   9.3242e-01
   9.3616e-01
   9.3986e-01
   9.4201e-01
   9.4145e-01
   9.4434e-01
   9.4252e-01
   9.4131e-01
   9.4249e-01
   9.4283e-01
   9.4395e-01
   9.4355e-01
   9.4690e-01
   9.4919e-01
   9.5255e-01
   9.5626e-01
   9.6282e-01
   9.6283e-01
   9.6672e-01
   9.6439e-01
   9.6887e-01
   9.7099e-01
   9.7529e-01
   9.8034e-01
   9.8302e-01
   9.8494e-01
   9.8828e-01
   9.8769e-01
   9.9130e-01
   9.9108e-01
   9.9422e-01
   9.9561e-01
   9.9737e-01
   1.0002e+00
   1.0034e+00
   1.0044e+00
   1.00329
   1.0];
d0 = 1.8e-10; % Initial guess for diffusion coefficient
fun = @(d) ypred(d, t1) - y_obs;
%******************************************************************************************************************************
%*****************************************************************************************************************************************
best_d = lsqnonlin(fun, d0) % Estimation of the best diffusion coefficient using lsqnonlin solver by calling the minimizing function and optimizing variable -d as arguments
predicted_y = ypred(best_d, t1); % Predicted At/Ainf curve corresponding to data fitting curve
% Plotting the lsqnonlin data fitting
figure(2)
plot(t1, y_obs, 'ro', t1, predicted_y, 'g-');
 
legend({'observed', 'predicted'})
ylabel('At/Ainf') 
xlabel('Time in seconds')
title('lsqnonlin data fitting')
%************
%Verification of whether the “lsqnonlin” solver accurately identified the local minimum to optimize the diffusion coefficient
pfun=@(d) norm(fun(d));
figure(1)
fplot(pfun,[0,5e-10])
hold on; plot(best_d*[1,1],ylim,'--rx'); 
title('Verification of whether minimum was correctly found')
hold off
%**********
% Diffusion model equation under function
function y_pred = ypred(d, t1)
    a=0.0011; % height of the tissue
    
    % The diffusion equation has been broken down for ease of tranformation
    % into code
    gama = 0.01005; 
    
    L2 = zeros(14,1);
    L3 = zeros(100,1);
    L4 = zeros(100,1);
    L5 = zeros(100,1);
    S= zeros(73,1);
    y_pred = zeros(73,1);
    
       
    L1 = ((8*gama)/((pi*(1-exp(-2*gama*a)))));
    format longE
    
   
    for t = t1(:).'
        for n=0:1:100
            L2(n+1) = exp((((2*n + 1)^2)*-d*pi*pi*t)/(4*a*a));
            L3(n+1) = (((-1)^n)*2*gama)+(((2*n+1)*pi)*exp(-2*gama*a))/(2*a);
            L4(n+1)= ((2*n)+1)*((4*gama*gama)+((((2*n)+1)*pi)/(2*a))^2);
            L5(n+1) = ((L2(n+1)*L3(n+1))/L4(n+1));
        end
        S((t/300) +1) = sum(L5);
        y_pred((t/300)+1)= 1 -(L1*S((t/300) +1)); % predicted data (before fitting)
    end
end

Fminbnd

% Non-linear data fitting to estimate the diffusion coefficent
% Solver = fminbnd
t1 = [0:300:21600]';
y_obs = [         
      0
   1.6216e-01
   2.9813e-01
   4.0805e-01
   4.9338e-01
   5.5928e-01
   6.0894e-01
   6.5506e-01
   6.8876e-01
   7.1773e-01
   7.4284e-01
   7.6433e-01
   7.8448e-01
   7.9912e-01
   8.1484e-01
   8.2950e-01
   8.3760e-01
   8.4707e-01
   8.5767e-01
   8.6299e-01
   8.6862e-01
   8.7433e-01
   8.7879e-01
   8.8736e-01
   8.9152e-01
   8.9903e-01
   9.0343e-01
   9.0984e-01
   9.1344e-01
   9.1615e-01
   9.2046e-01
   9.2313e-01
   9.2640e-01
   9.2992e-01
   9.3134e-01
   9.3242e-01
   9.3616e-01
   9.3986e-01
   9.4201e-01
   9.4145e-01
   9.4434e-01
   9.4252e-01
   9.4131e-01
   9.4249e-01
   9.4283e-01
   9.4395e-01
   9.4355e-01
   9.4690e-01
   9.4919e-01
   9.5255e-01
   9.5626e-01
   9.6282e-01
   9.6283e-01
   9.6672e-01
   9.6439e-01
   9.6887e-01
   9.7099e-01
   9.7529e-01
   9.8034e-01
   9.8302e-01
   9.8494e-01
   9.8828e-01
   9.8769e-01
   9.9130e-01
   9.9108e-01
   9.9422e-01
   9.9561e-01
   9.9737e-01
   1.0002e+00
   1.0034e+00
   1.0044e+00
   1.00329
   1.0];
%************
% Guessing the initial assumption d0 by finding the minimum using min function
%Implementing fminbnd 
pfun = @(d) norm( ypred(d, t1) - y_obs); %% Function whose sum of squares is minimized. 
dsamps=linspace(0,5e-10,50); %Generating a vector to plot above fubction and  identify the suitable initial guess
[~,imin]=min(  arrayfun(pfun,dsamps) ); % identification of the local minimum
[best_d,fval,exitflag,output]=fminbnd(pfun,dsamps(imin-1), dsamps(imin+1),...
                              optimset('TolX',1e-14)) % fminbnd solver to perform data fitting using the initial guess corresponding to the local minimum
exitflag,
best_d,
%****************
%Plot
predicted_y = ypred(best_d, t1);
figure(2)
plot(t1, y_obs, 'ro', t1, predicted_y, 'b-');
ylabel('At/Ainf') 
xlabel('Time in seconds')
legend({'observed', 'predicted'})
title('fminbnd solver data fitting')
%*************************************************************************
%Verification of whether the "fminbnd" solver accurately identified the local minimum to optimize the diffusion coefficient
figure(1)
fplot(pfun,[0,4e-10])
hold on; plot(best_d*[1,1],ylim,'--rx'); 
title('Verification of whether minimum was correctly found by fminbnd solver')
hold off
%*************************************************************************
% Diffusion model equation under function
function y_pred = ypred(d, t1)
    a=0.0011; % height of the tissue
    
    % The diffusion equation has been broken down for ease of tranformation
    % into code
    gama = 0.01005;
    
    L2 = zeros(14,1);
    L3 = zeros(100,1);
    L4 = zeros(100,1);
    L5 = zeros(100,1);
    S= zeros(73,1);
    y_pred = zeros(73,1);
    
    
    
    L1 = ((8*gama)/((pi*(1-exp(-2*gama*a)))));
    format longE
    
  
    for t = t1(:).'
        for n=0:1:100
            L2(n+1) = exp((((2*n + 1)^2)*-d*pi*pi*t)/(4*a*a));
            L3(n+1) = (((-1)^n)*2*gama)+(((2*n+1)*pi)*exp(-2*gama*a))/(2*a);
            L4(n+1)= ((2*n)+1)*((4*gama*gama)+((((2*n)+1)*pi)/(2*a))^2);
            L5(n+1) = ((L2(n+1)*L3(n+1))/L4(n+1));
        end
        S((t/300) +1) = sum(L5);
        y_pred((t/300)+1)= 1 -(L1*S((t/300) +1)); % predicted data (before fitting)
    end
end

Anand Ra 2021-7-2

Oh okay, got it. The percentage difference between the estimates is concerning. Looking at the curve, it looks like it’s un constrained in the end and probably the predicted curve is not even reaching the equilibrium. Is it possible to truncate the prediction by giving upper bounds to achieve that fit? I can try that actually but wanted to understand from you if the solvers can do that or misinterpret it.