How to curve fit an equation that gets values from another equation?

I have 38 data points for x and t (input variables) and 38 data points for y (output variable). In equation 1, a and b are coefficients to be found. z is the output variable which iterates and gets input from its previous value. That means when i=1, z_2 gets value from z_1. when i=2, z_3 gets value from z_2, and so on.

We can take z_0 = 0.5

.................equation (1)

solving equation 1 should give me 38 points for z, which will depend on the x and t.

....................equation (2)

In equation 2, x is the input data, y(fit) is the output and z comes from equation 1. c,d and e are coefficient to be found.

Now, I want to fit equation 2 to y vs x curve (which I have from experiments) and find the values of the constants a,b,c,d and e.

I have previously used lsqcurvefit to fit a single equation. But in this case there are 2 equations and the second one gets values from the first one.

Any help with this will be much appreciated. Thanks in advance.

[Rik]

You can write a custom function that does all calculations you need. What have you tried?

If you need help, please attach the relevant data. It is also not entirely clear to me which exact parameters need to be fitted and what data would be the input and output.

[Alex Sha]

It is right, give out the data of x, t and y, this will let others to help you more convenient。

[Ases Akas Mishra]

Thanks for your interest. I have attached an excel file which has x,t and y.

More clarification-

x and t are input data and y is the output data. In equation 1, we input x and t to get z. And in equation 2, we input x and the value of z from equation 1 to get y(fit). In the end, the goal is to fit y(fit) vs x curve to y vs x curve(from excel file). As equation 2 depends on z(from equation 1), we need to optimize the model in such a way that we obtain unique values of fit parameters (a,b,c,d,e) which satisfies both equation and the resulting curve fits y vs x data.

Till now I am manually supplying guessed values of a,b,c,d,e to obtain a curve. But in reality we need to get the values of these constants using some fitting function like lsqcurvefit.

I hope this is clear. Please let me know if you need more clarification.

[Rik]

First you need to load the data to Matlab variables.

What code are you using to calculate the curve? It will probably be fairly easy to adapt it.

[Ases Akas Mishra]

Hi, I am new to MATLAB central and I dont know how to load the data. I have pasted the data below. Let me know if you need it in some other format.

x = 1, 1.44, 2.07, 2.98, 4.28, 6.16, 8.86, 12.7, 18.3, 26.4, 37.9, 54.6, 78.5, 113, 162, 234, 336, 483, 695, 1000, 1000, 695, 483, 336, 234, 162, 113, 78.5, 54.6, 37.9, 26.4, 18.3, 12.7, 8.86, 6.16, 4.28, 2.98, 2.07

t = 1800, 3600, 5400, 7200, 9000, 10800, 12600, 14400, 16200, 18000, 19800, 21600, 23400, 25200, 27000, 28800, 30600, 32400, 34200, 36000, 37800, 39600, 41400, 43200, 45000, 46800, 48600, 50400, 52200, 54000, 55800, 57600, 59400, 61200, 63000, 64800, 66600, 68400

y = 1.419, 1.3369, 1.1017, 1.4513, 2.2703, 3.1121, 4.1855, 4.2835, 4.6603, 4.8933, 5.287, 5.877, 6.4844, 7.2159, 7.9723, 12.472, 15.202, 21.45, 24.348, 27.307, 26.128, 21.671, 18.912, 16.709, 14.279, 12.437, 10.866,9.5601, 8.5242, 7.733, 7.0755, 6.5612, 6.1244, 5.531, 4.8106, 3.2184, 2.1077, 1.6007

And here is my code-

clc
clear all
close all

x= xlsread('data.xlsx','0.5','B4:B41');
y= xlsread('data.xlsx','0.5','C4:C41');
t = xlsread('data.xlsx','0.5','H4:H41');

for i = 1 : length(x)
    
   % to find z
   
    a = 100;                         % I have used some random value for a and b. But they should be obtained from the fit.
    b = 0.05;
    
    A(i,:) = a/((b.*x(i)) + a);     % parts of equation 1
    e(i,:) = -t(i) / (a.*x(i)+b);
  
    z(1,1) = 0.5;                    %initialize z
    z(i+1,:) = A(i) + (z(i) - A(i)) * exp(e(i));   %equation 1
    z1{i,:} = z(i+1,:)                %to remove the extra term(z=0.5)
    z_final = cell2mat(z1);            %converting to double
    
    % equation 2
    
    c = 0.5;                    %Again guessing values for c,d,e
    d = 2.2;
    e = 0.35;
    
    y_fit(i,:) = (c + d.*z_final(i)).*x(i)^e   %equation 2
   
end

Answer 1

[Rik]

There is probably something wrong with your loop, as you're using e both in the first part as an array, and in the second part. To avoid errors I renamed the one in the first part to E. Since you're only ever accessing e(i) anyway, there is also no need for indexing (same for A).

This is all assuming your function actually works as expected. Tweak your function as needed.

x = [1, 1.44, 2.07, 2.98, 4.28, 6.16, 8.86, 12.7, 18.3, 26.4, 37.9, 54.6, 78.5, 113, 162, 234, 336, 483, 695, 1000, 1000, 695, 483, 336, 234, 162, 113, 78.5, 54.6, 37.9, 26.4, 18.3, 12.7, 8.86, 6.16, 4.28, 2.98, 2.07];
t = [1800, 3600, 5400, 7200, 9000, 10800, 12600, 14400, 16200, 18000, 19800, 21600, 23400, 25200, 27000, 28800, 30600, 32400, 34200, 36000, 37800, 39600, 41400, 43200, 45000, 46800, 48600, 50400, 52200, 54000, 55800, 57600, 59400, 61200, 63000, 64800, 66600, 68400];
y = [1.419, 1.3369, 1.1017, 1.4513, 2.2703, 3.1121, 4.1855, 4.2835, 4.6603, 4.8933, 5.287, 5.877, 6.4844, 7.2159, 7.9723, 12.472, 15.202, 21.45, 24.348, 27.307, 26.128, 21.671, 18.912, 16.709, 14.279, 12.437, 10.866,9.5601, 8.5242, 7.733, 7.0755, 6.5612, 6.1244, 5.531, 4.8106, 3.2184, 2.1077, 1.6007];

%initial estimates
a = 100;
b = 0.05;
c = 0.5;
d = 2.2;
e = 0.35;

lsqcurvefit(@(abcde,x)find_y_fit(x,t,abcde(1),abcde(2),abcde(3),abcde(4),abcde(5)),...
    [a b c d e],x,y)

Solver stopped prematurely.

lsqcurvefit stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 5.000000e+02.

ans = 1×5

83.0941 -0.0005 -337.2435 339.7435 0.2726

function y_fit=find_y_fit(x,t,a,b,c,d,e)
%it is possible to optimize/clean this a lot futher, but let's first get it working
y_fit=zeros(size(x));
z=zeros(size(x));
z(1,1) = 0.5; %initialize z
for n = 1 : numel(x)
   % to find z
    A = a/((b.*x(n)) + a);     % parts of equation 1
    E = -t(n) / (a.*x(n)+b);
  
    z(n+1,:) = A + (z(n) - A) * exp(E);   %equation 1
    z1{n,:} = z(n+1,:);                %to remove the extra term(z=0.5)
    z_final = cell2mat(z1);            %converting to double
    
    % equation 2
    y_fit(n) = (c + d.*z_final(n)).*x(n)^e ;  %equation 2 
end
end

[Ases Akas Mishra]

Hi Rik

Thanks for the solution.

I have the following issues and questions.

1. When I try to plot y vs x and y_fit vs x in the same plot I get the attached figure. I would expect the final fit line to appear in the plot and not all the iterations.
2. As you can see, the fit is not very good. This is because the initial values of a,b,c,d and e are not correct. Theoretically there exists an unique combination of a,b,c,d,e that gives perfect fit. How can I change the code so that MATLAB finds the initial values on it's own without e specifying it? I found something called 'Latin-hypercube-sampling' that does this operation. The function is 'lhsdesign' in MATLAB. Any suggestions on this?

I=imread('y_vs_x_rik.jpg');
imshow(I)

[Ases Akas Mishra]

Follow up request- I would also like to collect the final values of z_final and y_fit (after optimization) so that I can use it for other purposes such as plotting using another software.

[Rik]

I have little experience with advanced curve fitting. My naive approach would be to use something like ndgrid to generate a few initial estimates for the 5 parameters. Note that this will blow up quickly, as you will have to run the fit for N^5 iterations for N points per parameter.

Since you didn't provide the code you used, I can't tell whether the points simply lack sorting, or if there is something else wrong. You have not confirmed my edits to your code are actually correct. You really should. Use the debugger to step through the code line by line.

Answer 2

[Mathieu NOE]

hello

my first job was to make your code a bit more efficient and streamlined

also notice you use both "e" for variable name and for a constant , avoid this as it can lead to errors and unpredictable results. I renamed the "e" from your equation 1 as ""ee" to avoid any conflict with the constant "e".

so bellow is basically your code a bit rewanped :

clc
clear all
close all
% x= xlsread('data.xlsx','0.5','B4:B41');
% y= xlsread('data.xlsx','0.5','C4:C41');
% t = xlsread('data.xlsx','0.5','H4:H41');

data= xlsread('Data.xlsx');
x = data(:,1);
t = data(:,2);
y = data(:,3);

z = 0.5;
a = 100;                         % I have used some random value for a and b. But they should be obtained from the fit.
b = 0.05;
c = 0.5;                    %Again guessing values for c,d,e
d = 2.2;
e = 0.35;
    
for i = 1 : length(x)
    A(i,:) = a/((b.*x(i)) + a);     % parts of equation 1
    ee(i,:) = -t(i) / (a.*x(i)+b);
    z(i+1,:) = A(i) + (z(i) - A(i)) * exp(ee(i));   %equation 1
end

    z_final = z(2:end);            %to remove the (first) extra term(z=0.5)
    % equation 2
    y_fit = (c + d.*z_final).*(x.^e);   %equation 2
    
plot(t,y,t,y_fit);

second step, your equations are now a function to minimize the distance (norm) vs the experimental data , and you get the optimal parameters in array opt_vector ( that contains your five constants [a b c d e] ):

here I used the fminsearch function, but you can use more advanced optmization function if you need to add some constraints on the parameters (I don't have the Optimization Toolbox);

results :

opt_vector = 12.3055 -0.0016 -0.3599 1.9294 0.3889

code :

clc
clear all
close all
% x= xlsread('data.xlsx','0.5','B4:B41');
% y= xlsread('data.xlsx','0.5','C4:C41');
% t = xlsread('data.xlsx','0.5','H4:H41');

data= xlsread('Data.xlsx');
x = data(:,1);
t = data(:,2);
y = data(:,3);

%% fit equation to data
% parameters
a = 100;                         % I have used some random value for a and b. But they should be obtained from the fit.
b = 0.05;
c = 0.5;                    %Again guessing values for c,d,e
d = 2.2;
e = 0.35;
params = [a b c d e];
    
opt_vector = fminsearch(@ (params) norm (y - myequation(params, data)) ,params);

y_fitted = myequation(opt_vector, data);

figure(1);  
plot(t,y,t,y_fitted);
legend('raw','fit');


function y_fit = myequation(params, data)

x = data(:,1);
t = data(:,2);
y = data(:,3);

% parameters
a = params(1);                        
b = params(2);
c = params(3);                   
d = params(4);
e = params(5);
% init z
z = 0.5;

    for i = 1 : length(x)

        A(i,:) = a/((b.*x(i)) + a);     % parts of equation 1
        ee(i,:) = -t(i) / (a.*x(i)+b);
        z(i+1,:) = A(i) + (z(i) - A(i)) * exp(ee(i));   %equation 1

    end

    z_final = z(2:end);            %to remove the (first) extra term(z=0.5)
    
    % equation 2
    y_fit = (c + d.*z_final).*(x.^e);   %equation 2
    
end

[Mathieu NOE]

Just for information, I added a line to force the parameters to remain positive (I guessed that from the fact that your initial values are all positive)

so in the function code , I added :

params = abs(params); % make sure they stay always positive

I recognize this is not the best method and you should probably better use aan optimizer that handles the constraints better than my trick;

new parameters :

opt_vector = 51.0698 0.0000 0.7150 0.5857 0.4359

seems like several parameters solutions lead to very similar fitted plots

new plot :

new code :

clc
clear all
close all
% x= xlsread('data.xlsx','0.5','B4:B41');
% y= xlsread('data.xlsx','0.5','C4:C41');
% t = xlsread('data.xlsx','0.5','H4:H41');

data= xlsread('Data.xlsx');
x = data(:,1);
t = data(:,2);
y = data(:,3);

%% fit equation to data
% parameters
a = 100;                         % I have used some random value for a and b. But they should be obtained from the fit.
b = 0.05;
c = 0.5;                    %Again guessing values for c,d,e
d = 2.2;
e = 0.35;
params = [a b c d e];
    
opt_vector = fminsearch(@ (params) norm (y - myequation(params, data)) ,params);

y_fitted = myequation(opt_vector, data);

figure(1);  
plot(t,y,t,y_fitted);
legend('raw','fit');


function y_fit = myequation(params, data)

x = data(:,1);
t = data(:,2);
y = data(:,3);

% parameters
params = abs(params); % make sure they stay always positive
a = params(1);                        
b = params(2);
c = params(3);                   
d = params(4);
e = params(5);
% init z
z = 0.5;

    for i = 1 : length(x)

        A(i,:) = a/((b.*x(i)) + a);     % parts of equation 1
        ee(i,:) = -t(i) / (a.*x(i)+b);
        z(i+1,:) = A(i) + (z(i) - A(i)) * exp(ee(i));   %equation 1

    end

    z_final = z(2:end);            %to remove the (first) extra term(z=0.5)
    
    % equation 2
    y_fit = (c + d.*z_final).*(x.^e);   %equation 2
    
end

[Ases Akas Mishra]

Hi Mathieu

Thanks for the solution.

I have the following issues and questions.

1. When I try to plot y vs x and y_fit vs x in the same plot I get the attached figure. I would expect the final fit line to appear in the plot and not all the iterations.
2. As you can see, the fit is not very good. This is because the initial values of a,b,c,d and e are not correct. Theoretically there exists an unique combination of a,b,c,d,e that gives perfect fit. How can I change the code so that MATLAB finds the initial values on it's own without e specifying it? I found something called 'Latin-hypercube-sampling' that does this operation. The function is 'lhsdesign' in MATLAB. Any suggestions on this?

I=imread('mathieu.jpg');
imshow(I)

[Ases Akas Mishra]

Follow up request- I would also like to collect the final values of z_final and y_fit (after optimization) so that I can use it for other purposes such as plotting using another software.

Answer 3

[Alex Sha]

How about the results below:

Residual Sum of Squares (RSS): 34.9122402573179
Root of Mean Square Error (RMSE): 0.958510910040288
Correlation Coef. (R): 0.991557637426582
R-Square: 0.983186548338985

Parameter          	Best Estimate
--------------------	-------------
a	4.63375351494044E-8
b	-1.56073777954123E-7
c	1.80838380873546
d	-1.35268042423128
e	0.387487908425502

[Ases Akas Mishra]

Hi ALex

This looks fine. Like I have mentioned before, the issue is to somehow guess the initial values of a,b,c,d,e such that the fit is perfect. There will be an unique combination of a,b,c,d,e when the fit will be perfect and then the 'calculated y' curve won't be symmetric (along vertical axis), like it is now. Please share your code. Thanks for your interest.