how to perform data fit like excel? and plot
Mostra commenti meno recenti
- I have observed array of data ( y_obs) and predicted data (y_pred)
- Predicted data is obtained from an equation
- How do I fit the observed data to the predicted data by minimizing the coefficient "d" in the equation? ( This is possible in excel, but I could not find a suitable method in matlab
Below is my code for steps 1 and 2:
% observed data
y_obs = [0.3 0.2 0.28 0.318 0.421 0.492 0.572 0.55 0.63 0.61 0.73 0.8 0.81 0.84 0.93 0.91]'; % If y_obs should equal to predicted, I can have more data. J us fo rthe code, I am providing limited observed data
t1 = [300:300:21600]';
a=0.0011;
gama = 0.01005;
d=0.000000000302;
n=1;
t=300;
L2 = zeros(14,1);
L3 = zeros(14,1);
L4 = zeros(14,1);
At = zeros(14,1);
t = 300;
n =1;
L1 = ((8*gama)/((pi*(1-exp(-2*gama*a)))));
format shortE
for t= 300:300:21600
for n=1:1:14
L2(n) = exp((((2*n + 1)^2)*-d*pi*pi*t)/(4*a*a));
L3(n) = (((-1)^n)*2*gama)+(((2*n+1)*pi)*exp(-2*gama*a))/(2*a);
L4(n)= ((2*n)+1)*((4*gama*gama)+((((2*n)+1)*pi)/(2*a))^2);
L5(n) = ((L2(n).*L3(n))/L4(n));
end
S(t/300) = sum(L5);
y_pred(t/300) = 1 -L1*S(t/300); % predicted data
end
2 Commenti
Walter Roberson
il 16 Giu 2021
L2 = zeros(14);
that should probably be
L2 = zeros(14,1);
like the other variables.
Anand Ra
il 16 Giu 2021
Risposta accettata
Walter Roberson
il 16 Giu 2021
Modificato: Walter Roberson
il 16 Giu 2021
format shortE
% observed data
y_obs = [0.3 0.2 0.28 0.318 0.421 0.492 0.572 0.55 0.63 0.61 0.73 0.8 0.81 0.84 0.93 0.91]'; % If y_obs should equal to predicted, I can have more data. J us fo rthe code, I am providing limited observed data
t1 = [300:300:21600]';
T1 = t1(1:length(y_obs)).';
a=0.0011;
gama = 0.01005;
d0 = 0.000000000302;
syms d
n=1;
t=300;
L2 = sym(zeros(14,1));
L3 = sym(zeros(14,1));
L4 = sym(zeros(14,1));
At = sym(zeros(14,1));
t = 300;
n =1;
L1 = ((8*gama)/((pi*(1-exp(-2*gama*a)))));
y_pred = sym(zeros(length(T1),1));
for t = T1
for n=1:1:14
L2(n) = exp((((2*n + 1)^2)*-d*pi*pi*t)/(4*a*a));
L3(n) = (((-1)^n)*2*gama)+(((2*n+1)*pi)*exp(-2*gama*a))/(2*a);
L4(n)= ((2*n)+1)*((4*gama*gama)+((((2*n)+1)*pi)/(2*a))^2);
L5(n) = ((L2(n).*L3(n))/L4(n));
end
S(t/300) = sum(L5);
y_pred(t/300) = 1 -L1*S(t/300); % predicted data
end
sse = expand(sum((y_pred - y_obs(:)).^2));
f = matlabFunction(sse)
opt1 = fmincon(f, d0)
opt2 = fminsearch(f, d0)
35 Commenti
Walter Roberson
il 16 Giu 2021
My test system got really slow suddenly, so I have not tested the fmincon / fminsearch steps yet.
Walter Roberson
il 16 Giu 2021
Please note that you are building 72 y_pred but your y_obs only has 16 entries. I coded above to only go as far in predictions as you have y_obs for.
Walter Roberson
il 16 Giu 2021
fmincon suggests d about 1e-11 and fminsearch suggests d about 3e-13 . But d = 0 is better than 1e-11, and the value at -3e-13 is a little better than that. -3.3508267261899e-13 checks out as a local minima when I plot.
Anand Ra
il 16 Giu 2021
Walter Roberson
il 16 Giu 2021
L5 = sym(zeros(14,1));
right after assignment to L4
Anand Ra
il 16 Giu 2021
Anand Ra
il 16 Giu 2021
Anand Ra
il 17 Giu 2021
Walter Roberson
il 17 Giu 2021
%assuming that opt is the optimal d value:
y_pred_n = double(subs(y_pred, d, opt));
plot(T1, y_obs, 'k*', T1, y_pred_n, 'b')
legend({'observed', 'modeled'})
Anand Ra
il 17 Giu 2021
Anand Ra
il 17 Giu 2021
format long g
% observed data
y_obs = [0.3 0.2 0.28 0.318 0.421 0.492 0.572 0.55 0.63 0.61 0.73 0.8 0.81 0.84 0.93 0.91]'; % If y_obs should equal to predicted, I can have more data. J us fo rthe code, I am providing limited observed data
t1 = [300:300:21600]';
T1 = t1(1:length(y_obs)).';
a=0.0011;
gama = 0.01005;
d0 = 0.000000000302;
syms d
n=1;
t=300;
L2 = sym(zeros(14,1));
L3 = sym(zeros(14,1));
L4 = sym(zeros(14,1));
At = sym(zeros(14,1));
t = 300;
n =1;
L1 = ((8*gama)/((pi*(1-exp(-2*gama*a)))));
y_pred = sym(zeros(length(T1),1));
for t = T1
for n=1:1:14
L2(n) = exp((((2*n + 1)^2)*-d*pi*pi*t)/(4*a*a));
L3(n) = (((-1)^n)*2*gama)+(((2*n+1)*pi)*exp(-2*gama*a))/(2*a);
L4(n)= ((2*n)+1)*((4*gama*gama)+((((2*n)+1)*pi)/(2*a))^2);
L5(n) = ((L2(n).*L3(n))/L4(n));
end
S(t/300) = sum(L5);
y_pred(t/300) = 1 -L1*S(t/300); % predicted data
end
sse = expand(sum((y_pred - y_obs(:)).^2));
f = matlabFunction(sse)
[opt_mc, err_mc] = fmincon(f, d0)
[opt_ms, err_ms] = fminsearch(f, d0)
y_pred_n = double(subs(y_pred, d, opt_mc));
plot(T1, y_obs, 'k*', T1, y_pred_n, 'b')
legend({'observed', 'modeled'})
title('fmincon predictions')
y_pred_n = double(subs(y_pred, d, opt_ms));
plot(T1, y_obs, 'k*', T1, y_pred_n, 'b')
legend({'observed', 'modeled'})
title('fminsearch predictions')
These are not straight lines, and the predicted minimum is not the same as d0.
The visual fit is not good in either case.
Walter Roberson
il 17 Giu 2021
If you continue on with
trypoints = -logspace(-14,-5, 10000);
predictions = f(trypoints);
[bestprediction, bestidx] = min(predictions)
bestd = trypoints(bestidx)
trypointsp = logspace(-14,-5, 10000);
predictions = f(trypointsp);
[bestpredictionp, bestidx] = min(predictions)
bestdp = trypointsp(bestidx)
then the bestd is still around -3.34774581070573e-13
Either your calculation is quite sensitive to numeric values (that is, if you were to proceed step by step numerically rather than creating a formula and substituting into that, then you might get a different result due to round-off happening much more often)... or else the model is not a good fit for the data.
Anand Ra
il 18 Giu 2021
Walter Roberson
il 18 Giu 2021
Constraints will not help much. All a constraint could help you do would be to avoid going below about -4e-13, which is the point at which it becomes NaN.
My tests show that at least for the data that you posted earlier, that the shape goes:
- almost everything negative gives NaN
- between roughly -4e-13 and -3.35e-13, the error calculation exists and decreases like a quadratic
- above roughly -3.35e-13 and 0, the value increases again like a quadratic
- above 0 the value keeps increasing
When I took the derivative of y_pred and solved for it being 0, the solution was the approximately -3.35 value.
If your assignment statements are correct, then the model is simply wrong for the data. Which implies you should check your assignment statements in detail.
Anand Ra
il 18 Giu 2021
format long g
% observed data
y_obs = [0.3 0.2 0.28 0.318 0.421 0.492 0.572 0.55 0.63 0.61 0.73 0.8 0.81 0.84 0.93 0.91]'; % If y_obs should equal to predicted, I can have more data. J us fo rthe code, I am providing limited observed data
t1 = [300:300:21600]';
T1 = t1(1:length(y_obs)).';
a=0.0011;
gama = 0.01005;
d0 = 0.000000000302;
syms d
n=1;
t=300;
L2 = sym(zeros(14,1));
L3 = sym(zeros(14,1));
L4 = sym(zeros(14,1));
At = sym(zeros(14,1));
t = 300;
n =1;
L1 = ((8*gama)/((pi*(1-exp(-2*gama*a)))));
y_pred = sym(zeros(length(T1),1));
for t = T1
for n=1:1:14
L2(n) = exp((((2*n + 1)^2)*-d*pi*pi*t)/(4*a*a));
L3(n) = (((-1)^n)*2*gama)+(((2*n+1)*pi)*exp(-2*gama*a))/(2*a);
L4(n)= ((2*n)+1)*((4*gama*gama)+((((2*n)+1)*pi)/(2*a))^2);
L5(n) = ((L2(n).*L3(n))/L4(n));
end
S(t/300) = sum(L5);
y_pred(t/300) = 1 -L1*S(t/300); % predicted data
end
flsq = matlabFunction(y_pred(:))
[opt_lsq, err_lsq] = lsqnonlin(flsq, d0)
y_pred_n = double(subs(y_pred, d, opt_lsq));
plot(T1, y_obs, 'k*', T1, y_pred_n, 'b')
legend({'observed', 'modeled'})
title('lsqnonlin predictions')
Walter Roberson
il 18 Giu 2021
I recommend that you take an approach like I showed in https://www.mathworks.com/matlabcentral/answers/857320-how-to-perform-data-fit-like-excel-and-plot#comment_1589825 where you construct a vector of trial d values to evaluate f at, and plot the result. You will be able to see visually that there isn't any hidden low-residue area, that you are not dealing with multiple minima and maxima. Not unless mabye you start going for much higher d values... I can't say that I plotted for (say) d = 317.
Anand Ra
il 19 Giu 2021
Your model cannot even get close if you confine to two observations :(
If you confine to one observation, then the fminsearch() version can get quite close for that one point, but the fmincon() version doesn't get close.
Taylor series doesn't help it seems.
format long g
% observed data
y_obs = [0.3 0.2 0.28 0.318 0.421 0.492 0.572 0.55 0.63 0.61 0.73 0.8 0.81 0.84 0.93 0.91]'; % If y_obs should equal to predicted, I can have more data. J us fo rthe code, I am providing limited observed data
t1 = [300:300:21600]';
T1 = t1(1:length(y_obs)).';
y_obsX = y_obs([1 3]);
T1X = T1([1 3]);
a=0.0011;
gama = 0.01005;
d0 = 0.000000000302;
syms d
n=1;
t=300;
L2 = sym(zeros(14,1));
L3 = sym(zeros(14,1));
L4 = sym(zeros(14,1));
At = sym(zeros(14,1));
t = 300;
n =1;
L1 = ((8*gama)/((pi*(1-exp(-2*gama*a)))));
y_pred = sym(zeros(length(T1X),1));
for tidx = 1 : length(T1X)
t = T1X(tidx);
for n=1:1:14
L2(n) = exp((((2*n + 1)^2)*-d*pi*pi*t)/(4*a*a));
L3(n) = (((-1)^n)*2*gama)+(((2*n+1)*pi)*exp(-2*gama*a))/(2*a);
L4(n)= ((2*n)+1)*((4*gama*gama)+((((2*n)+1)*pi)/(2*a))^2);
L5(n) = ((L2(n).*L3(n))/L4(n));
end
S(tidx) = sum(L5);
y_pred(tidx) = 1 - L1*S(tidx); % predicted data
end
sse = expand(sum((y_pred - y_obsX(:)).^2));
f = matlabFunction(sse)
[opt_mc, err_mc] = fmincon(f, d0)
[opt_ms, err_ms] = fminsearch(f, d0)
y_pred_n = double(subs(y_pred, d, opt_mc));
plot(T1, y_obs, 'k*', T1X, y_pred_n, 'r-v')
legend({'observed', 'modeled'})
title('fmincon predictions')
y_pred_n = double(subs(y_pred, d, opt_ms));
plot(T1, y_obs, 'k*', T1X, y_pred_n, 'r-v')
legend({'observed', 'modeled'})
title('fminsearch predictions')
best_d = vpasolve(diff(sse, d))
f(double(best_d))
FT = taylor(sse, d, 0, 'order', 10)
trial_sol = vpasolve(FT)
f(double(trial_sol))
Walter Roberson
il 20 Giu 2021
Above d=1e-9, your function pretty much predicts 1.0 for all values
There just isn't anywhere left to go. Your y_pred do not model your data.
Anand Ra
il 20 Giu 2021
Anand Ra
il 20 Giu 2021
Walter Roberson
il 21 Giu 2021
You can perform fitting with numeric data.
t1 = [0:300:21600]';
y_obs = [
0
1.6216e-01
2.9813e-01
4.0805e-01
4.9338e-01
etc];
d0 = 0.000000000302;
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, '-');
function y_pred = ypred(d, t1)
%HERE
end
at the point marked HERE insert your code that takes a numeric d and numeric t1 and calculates y_pred -- basically just moving some of your code around, such as moving a and gama inside of the new function, and otherwise mostly doing the same thing
Anand Ra
il 21 Giu 2021
Anand Ra
il 21 Giu 2021
I plotted the equation predicted vs solver predicted. Just overlaps. No minimation is perfpormed by the solver
solver predicted
solver predicted
Walter Roberson
il 21 Giu 2021
Modificato: Walter Roberson
il 21 Giu 2021
You overwrote d and t1 in the function and you hardcoded t1 in the function.
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 = 1e-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
Walter Roberson
il 21 Giu 2021
But still notice that the fitting is giving the same d result as input. The fitting can only progress down to about 1e-8 with the default fitting options.
Anand Ra
il 21 Giu 2021
Anand Ra
il 25 Giu 2021
Hello,
All of a sudden, the curvesa are fitting. I have no idea what I did.
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 = 3.0e-10;
fun = @(d) ypred(d, t1) - y_obs;
best_d = lsqnonlin(fun, d0)
predicted_y = ypred(best_d, t1);
plot(t1, y_pred, '*', 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);
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
il 25 Giu 2021
Anand Ra
il 26 Giu 2021
Più risposte (0)
Categorie
Scopri di più su Spectral Measurements in Centro assistenza e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
