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

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

[Matt J]

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]

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]

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.

[Anand Ra]

Okay, do you mean, I do perform the grid search over the y_obs or y_pred like below? And then pass into the function? (I am getting an error) or do that seperately, find d0 and then pass it into lsqnonlin?

[X Y] = ngrid(y_pred,t1(:)')

[d] = ngrid(0:1e-15:5e-10)

[Anand Ra]

OR is it easier to use fminsearch or lscurvefit to do the above? Thanks!

[John D'Errico]

Better to use fminbnd than either fminsearch or lsqcurvefit. fminbnd will be far better behaved.

[Anand Ra]

Sorry, I am having trouble understanding. Should I use fminsearch to a find the initial guess and integrate that to the lsqnonlin code? Or you are asking me to totally replace lsqnonlin. I have the first scenario coded below. Not sure if it makes sense, as its spitting out the same "d" everytime, irresptive of my limits.

function y_pred = best_d(d)

t1 = [0:300:21600]'; 
L1 = ((8*gama)/((pi*(1-exp(-2*gama*a)))));
    format longE

a=0.0011;
    gama = 0.01005;
    
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

d = fminbnd(@best_d,2e-10, 4e-10)

[Matt J]

For example, using your originally posted data set, I could do,

load Data

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,

exitflag = 
     1

best_d,

best_d = 
     1.626778798935391e-10

[Anand Ra]

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

[Matt J]

That definitely helps with the initial guess. However, wondering why I am getting even better fit if I give d0 = 3.2e-10

I'm definitely not seeing that with the data set you've given us.

[Anand Ra]

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

[Matt J]

The global minimum was found. Constraining the fit would just make it worse.

[Anand Ra]

Got it, thanks!! . So, then, is this is the best fit that can be obtained?

[Matt J]

With the model you are using, yes.

[Anand Ra]

Alright, Thank you Sir. I will try with different data sets to verify the fit.

[Anand Ra]

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

[Matt J]

It does look a bit different , but it isn't much of a surprise that a different model gives a different fit.

[Anand Ra]

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.

[Matt J]

There are tools that can give you better control over the shape of the fit, for example,

https://www.mathworks.com/matlabcentral/fileexchange/24443-slm-shape-language-modeling