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

조회 수: 20 (최근 30일)
Ases Akas Mishra
Ases Akas Mishra 2022년 2월 16일
댓글: Rik 2022년 2월 20일
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.
  댓글 수: 5
Rik
Rik 2022년 2월 17일
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
Ases Akas Mishra 2022년 2월 17일
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

댓글을 달려면 로그인하십시오.

답변 (3개)

Rik
Rik 2022년 2월 18일
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
  댓글 수: 3
Ases Akas Mishra
Ases Akas Mishra 2022년 2월 20일
편집: Ases Akas Mishra 2022년 2월 20일
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
Rik 2022년 2월 20일
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.

댓글을 달려면 로그인하십시오.


Mathieu NOE
Mathieu NOE 2022년 2월 18일
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
  댓글 수: 3
Ases Akas Mishra
Ases Akas Mishra 2022년 2월 20일
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
Ases Akas Mishra 2022년 2월 20일
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.

댓글을 달려면 로그인하십시오.


Alex Sha
Alex Sha 2022년 2월 18일
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
  댓글 수: 1
Ases Akas Mishra
Ases Akas Mishra 2022년 2월 20일
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.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Least Squares에 대해 자세히 알아보기

태그

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by