Curve fitting and convergence to estimate two coefficients
Show older comments
Hello
Question Objective:
I am looking to employ curve fiting and convergence to estimate two coefficients in my model equation. I have estimated one coefficient by lscurvefit, fminsearch, but wondering if there is any solver in Matlab to estimate two coefficients in the approach ( details below) that I prefer.
Model Equation: ( See the below code for reference)
- Model equation is constructed using the terms L1, L2, L3, L4, L5, S and At_Af
- Final form of the model:
At_Af(t+1) = Af*(1 - L1*S((t)+1)
Coefficients to be estimated:
- Two coefficients - d and Af
Apporach that I wanted to employ:
- Consider the base line data( curve) array for fitting
- Evaluate the best fit by stepping 'n' in the equation from 0 to #( how many ever iterations) to obtain the best fit while trying to determine the two coefficients d and Af
- Can consider d= 1*10^-10 and Af = 0.100826 as the initial guess for the regression.
Data set - attachment:
The initial data set ( baseline curve data) = At_Af (LHS of the model equation) is attached to this query in a spreadsheet
The code of the model equation is below:
L=0.00075; % height of the tissue
% The diffusion equation has been broken down for ease of tranformation
% into code
gama = 0.000167;
L2 = zeros(14,1);
L3 = zeros(100,1);
L4 = zeros(100,1);
L5 = zeros(100,1);
S= zeros(73,1);
At_Af = zeros(73,1);
t1 = 0:1:3000
d= 1*10^-10;
L1 = ((8*gama)/((pi*(1-exp(-2*gama*L)))));
format longE
t_min = t1./60
for t = t1(:).'
for n=0:1:50
L2(n+1) = exp((((2*n + 1)^2)*-d*pi*pi*t)/(4*L*L));
L3(n+1) = (((-1)^n)*2*gama)+(((2*n+1)*pi)*exp(-2*gama*L))/(2*L);
L4(n+1)= ((2*n)+1)*((4*gama*gama)+((((2*n)+1)*pi)/(2*L))^2);
L5(n+1) = ((L2(n+1)*L3(n+1))/L4(n+1));
end
S((t) +1) = sum(L5);
At_Af(t+1) = Af*(1 - L1*S((t)+1));
end
Any help would be greatly appreciated. Thanks in advance!
Accepted Answer
William Rose
on 22 Sep 2021
0 votes
I was bothered by the poor performance of fmincon() on this problem. Even when we tried starting at multiple locations, it did not progress reliably o the true minimum. I figured out why: the values of the parameter vector, x=[d,Af], differ by 10 orders of magnitude (d~=1e-10, Af~=1). This large difference in scale of the parameters is a problem for fmincon(). fmincon() works much better if we scale the elements of the vector x to have the same order of magnitude. Threfore I adapted the code so that we specify d*1e10 in the main program. d*1e10 has magnitude~=1, like Af. This is the value passed to fmincon(). We divide this number by 1e10 inside myModel(), to get it back to its "correct" value. With this adjustment, we don't need the "PolyStart" version of the program. You can use fitData.m, which makes only one initial guess. You will see that for a wide rage of initial guesses, it fiinds the same best fit result - which is good. And the best fit it find is a lot better than the "best fit" we found before: the rms error is 0.00000.
I changed the allowable range for Af to [0.5,1.5], compared to the range [0,1] which you specified. I did this because myModelTest.m shows that values of Af in the range [0.5 to 1.5] are reasonable.
The modified versions of fitData.m, fitDataPolyStart.m, and myModel.m, are attached.
Console output and graphical output below, when the initial guess is d0*1e10=.5, Af0=.5:
>> fitData
Best fit at d=1.000e-10, Af=1.0000. R.m.s.error=0.00000.
fitDataPolyStart produces identical reults, after trying 30 different starting guesses.
14 Comments
Anand Ra
on 22 Sep 2021
William Rose
on 22 Sep 2021
There was a line in myModel.m which was
for i=0:3000
I have changed it to
for i=0:N-1
I saved the new file (with 241 points) as AtAfExpt2.txt, and I changes the file name to AtAfExpt2.txt in fitData.m. With these changes, fitData.m runs without error.
But the fit is not very good:
>> fitData
Best fit at d=7.155e-11, Af=0.5227. R.m.s.error=0.01078.
The graph of experiment versus simulation above, and our understanding of d and Af, which we gained from running myModelTest.m, make it clear that we need a higher value of d and a smaller value of Af to improve the fit. This reminds me that I changed the bounds for Af to [.5,1.5]. I will therefore change the bounds for Af to [0,1.5] in fitData.m. With this change, fitData.m finds a much better fit - see below. The best fit value of d is higher now, and the value of Af is smaller. Those are the changes we predicted.
>> fitData
Best fit at d=5.047e-09, Af=0.0847. R.m.s.error=0.00176.
The updated versions of fitData.m and myModel.m are attached.
Anand Ra
on 22 Sep 2021
William Rose
on 22 Sep 2021
@Anand Rathnam, you're welcome.
Anand Ra
on 1 Oct 2021
Anand Ra
on 1 Oct 2021
William Rose
on 2 Oct 2021
@Anand Rathnam, You posted a data file and 3 scripts with the same names as files we echanged previously. Are all four files altered from the earlier files with the same names? If so, how? It takes me extra time to figure this out.
I have renamed the files you sent me. I have named the latest data file AtAfExpt3.txt on my computer, since we already had AtAfExpt2.txt. I renamed your latest main program fitdata2.m on my computer, to avoid confusion with the earlier version, which works. Renaming function files (for example by adding a number to the filename) can cause challenges, because if you give the file a new name, you must rename the function itself in line 1, and you must change the call to the function, in any script that calls it. Therefore I have not renamed the function files myModel.m and sumsqerr.m, but I have saved the previous versions, which worked, with a filename reflecting the last date edited.
William Rose
on 2 Oct 2021
@Anand Rathnam, I am looking at your new main program which I have rnamed fitData2.m.
You decreased the minimum value for d*1e10, the first filled parameter, from 0.001 to 0.0001. Why?
The range you specified for the second parameter* to be fitted is now 0 to 0.005, with an initial guess of 0. I said in my ealrier comments that fmincon() wirks better if the parameters to be fitted are of the same order of magnitude. Your new range does not satisfy this suggestion since the max value is .005, whereas the max value for d*1e10 is 1000, and the best fit value on earlier runs is 1.
*The range and the initial guess for the second variable to be fitted have been altered, but none of the comments have been altered. The comments in fitData2.m still refer to the second parameter as Af. Do you actually mean for the second fitteted variable to be "L" now, rather than Af? If so, then fix the comments in the code, all the way through.
The comments in myModel.m also refer to Af, but I think you mean L. If that is true, then update the comments for myModel.m also.
You modified a line and added a line to myMOdel.m to show progress of minimization. The printed line refer to d and Af. Do you mean d and L?
The new program runs to completion without errors, and I can confirmthat the fitted value of the socndparameter is 0, which equals the initial guess adnd it at the low end of the range. I can alsos confirm that the fit doen't look too good. Do you know that a fitted value of 0 for L is wrong? Maybe that is in fact the best possible value in the range.
William Rose
on 5 Oct 2021
I understand that:
- You have set Af to unity. SInce Af acted as a scaling fctor, setting it to unity means that it no longer appears in the model equations.
- You want to fit d and L instead of d and Af.
- You have made good changes to the code in order to fit d and Af.
- The code runs to completion without error.
- The best-fit value of d is 0.001, which is at the low extreme of the allowable range.
You said in your Oct 3 comment that, in the earlier version of the code, L=0.0023. That is not correct, at least not for the version you posted on Sept 22, which is the version I used as my base for fitting. In the Sept 22 versoin of myModel.m, L=0.00075. The allowed range for L, in th Oct 3 version of fitData.m, is 0.001<=L<=5000. L=0.001 is the closest allowed value to the original value.
In the Oct 3 which you posted, code, gama=5678. In the Sept 22 code, gama=0.000167. This change in gama by a factor of more than 1e7 means that the best-fit values are likely to change.
By the way, your initial value for L in fitdata.m (Oct 3 versioin) is L=0. This is outside the allowed range. I highly recommend that the initial value should be in the interior of the allowed range.
William Rose
on 5 Oct 2021
The data in AtAfExpt3.txt, which you posted on Oct 3, has a jump from t=0 to t=1,2,3 (row 1 to rows 2,3,4...) Is the value at t=0 (row 1) correct?
When I run your code with the new experimental data file, AtAfExpt3.txt, I get the following console output and plot:
Best fit at d=1.313e-08, L=0.0010, R.m.s.error=0.06173.
I get a better fit (i.e. lower RMSE) in Excel, by setting d=1e-10, and adjusting L to make the error small. When L=1.03e-4, RMSE=0.029, which is a lot better than the fit above. This is with gama=5678, since that is the value of gama in myModel.m, which you posted on Oct 3. See Excel plot below.
If I use the Solver function in Excel, I can do somewhat better. See Excel plot below, in which the RMSE is 0.0243, which is better. But the values of d and L are outside the range you like. I think the unexpected best-fit values for d and L are due to gama=5678.
Now back to Matlab: Let's set gama=0.000167 in myModel.m. That is the value gama had on Sept 22. Then run fitData.m, with 1e-5<=L<=100, and initial guess L0=0.01. fitData.m finds the following solution:
Best fit at d=4.250e-10, L=2.286e-04, R.m.s.error=0.02444.
Anand Ra
on 12 Oct 2021
William Rose
on 13 Oct 2021
You're welcome, Anand. You said "Very interesting to see how excel out performs matlab solver. Did you attempt fitting it manually in excel."
I would not say excel outperforms Matlab solver, because they are very different, and they may be compared with various metrics giving different results.
I did try doing the fit manually before remembering thatExcel has an optimization capability.
More Answers (5)
William Rose
on 20 Sep 2021
fmincon() is a multidimensional minimization routine. By multidimensional I mean it adjusts multiple unknown parameters until a specified function (such as sum of squared errors) is minimized. fmincon() is complicated, but the help page for it is good, and has a bunch of examples. You will want to write a function that implements your diffusion model. You will pass to it the vector x=[d,Af]. It will return the vector AtAf.
function AtAf = myModel(x);
%MYMODEL Function to simulate diffusion
%Input: x=[d,Af];
%Output: AtAf=vector with same length as the experimental data vector
N=99; %length of experimental data vector; replace 99 with correct value
d=x(1);
Af=x(2);
AtAf=d*ones(N,1)+Af; %replace this line with code for the diffusion model
end
Save the code above as file myModel.m.
Write a function to compute the error between the model and the experimental data:
function sse=sumsqerr(x)
%SUMSQERR Function to compute sum squared error between model and experiment
%Input: x=[d,Af];
%Output: sse=scalar
AtAfExp=load('MyExptData.txt'); %read experimental data from text file
AtAfSim=myModel(x);
sse=sum((AtAfExp-AtAfSim).^2); %sum squared error
end
Save code above as file sumsqerr.m.
Write the main program. It calls fmincon(). You pass to fmincon() the initial guess x0=[d0;Af0], and the name of the function to be minimized (sumsqerr), and the allowed bounds for d and Af.
%fitData.m
x0=[1E-10,0.100826]; %initial guess
lb=[0,0]; %lower bound for [x,Af], replace with appropriate values
ub=[1,1]; %upper bound for [x,Af]; replace with appropriate values
x=fmincon(sumsqerr,x0,[],[],[],[],lb,ub);
disp(x) %display value of x=[d,Af] that minimizes the error
Save code above as fitData.m. Run fitData.m.
That gives you the general idea. Obviousy you need to make adjustments in the code above.
William Rose
on 21 Sep 2021
0 votes
Your code was a very good start. There were two fatal errors. ONe is that when you call fmincon(), you must pass a handle to the function, so you write fmincon(x0,@sumsqerr,[],[],[],[],lb,ub). I forgot to include the "@" in my example. The other problem is that function myMOdel() did not assign a value to AtAf. It assigned a value to At_Af. So I fied that. I also simplified your code and I replaced a for n=0:50 loop with vectorized statements to assign values to L2, L3, L4, L5. I checked that my vectorized versions were equivalent to the for loop values before deleting the foor loop code. Another change I made is that I read in the value for AtAfExp from a text file (see attached text file) in the main program, rather than defining it on every call to sumsqerr(). I declare AtAfExp to be global, in the main program and in sumsqerr(), so that its value will be known inside sumsqerr(). This approach may be faster, and more importantly it allows me to plot the eperimental and best fit values of AtAf one the "best fit" is found.
Now let's consider the results. The plot of the eperimental data and the "best fit' is below. It does not look good. I also show y=sqrt(t)/68, which is much better fit than the "best fit" simulation.
I think the "best fit" is not a good fit because fmincon() has gotten stuck in a local minimum that is not a global minimum. This is not uncomon in multidimensinal fitting with fmincon() and similar routines.
What can we do about it? We can manually try different values for the initial guess, x0. Or we can write code that automatically tries different starting points and finds the best fit from each starting point, and then picks the overall best fit. This is the approach I always take when doing multidimensional minimization.
William Rose
on 21 Sep 2021
Here is a new main program fitDataPolyStart.m, that tries multiple starting points, to reduce the chance of finding a solution that is not a global minimum. Use this instead of fitData.m. The other scripts and functions which I posted have not been changed.
Recal that with fitData, we obtained the "best fit"
d=4.1e-10, Af=0.560, RMSE=0.1093.
My first attempt with the PolyStart code was to try initial guessess that were spread equally across the allowed ranges which you specified. You specified 1e-15<=d<=1e-5, and 0<=Af<=1. Thefore I specified
d0=[1e-15,1e-13,1e-11,1e-9,1e-7,1e-5]; %initial guesses for d
Af0=[0,.2,.4,.6,.8,1]; %initial guesses for Af
This makes 6x6=36 initial guesses. The best of all 36 fits was
d=3.5e-10, Af=0.599, RMSE=0.0878.
That is better but still not very good. Therefore I refined by grid of starting points as follows:
d0=[1e-11,2e-11,5e-11,1e-10,2e-10,5e-10,1e-9]; %initial guesses for d
Af0=[.6,.7,.8,.9,1]; %initial guesses for Af
This makes 7x5=35 initial guesses. The best of all 35 fits was
d=1.8e-10, Af=0.775, RMSE=0.0264.
This is better. Refine the initial guesses again. I think we need a smaller d value and a bigger Af value to improve the fit.
d0=[1.5e-11,2e-11,3e-11,5e-11,7e-11,1e-10,1.5e-10,2e-10]; %initial guesses for d
Af0=[.75,.8,.85,.9,.95,1]; %initial guesses for Af
This best fit of these 56 initial guesses is at
d=1.50e-10, Af=0.850, RMSE=0.0139.
Notice that the RMSE improved by about a factor of 2.
Further guessing produces
d0=[1.0e-10,1.1e-10,1.2e-10]; %initial guesses for d
Af0=[.95,.96,.97]; %initial guesses for Af
which gives a very nice fit:
d=1.10e-10, Af=0.9600, RMSE=0.0026.
I have had good results with fmincon() for various model fitting problems in the past. However, I don't think fmincon() is doing a very good job on this fitting problem. IIt seems to get stuck, or it stops tryng to improve, too soon. fmincon() has many options. Perhpas the performance would improve with different choices for certain options.
William Rose
on 21 Sep 2021
0 votes
Here is fitDataPolyStart.m.
William Rose
on 21 Sep 2021
0 votes
I also wrote the attached script: myModelTest.m. This script calls test myModel to make sure it gives reasonable output, and to see how changes in d and Af affect the output. It does not do a search or a minimization. It just sets the values for d and Af and computes the corresponding AfAt(t). It does it for five different combinations of d and Af. I like to choose a "central" pair of [d,Af]. Then I do higher and lower d, with the same Af; and higher and lower Af, with the same d. That lets me see what Af does to the fit, and what d does to the fit.
The script also reads the experimental data from a file, and it computes the RMS Error of each of the five model outputs, reative to the experimental data.
I learned from this that d is like a rate constant for exponential rise or decay. When d is big, AfAt(t) quickly approaches its asymptotic value. When d is small, AfAt(t) takes a long time to approach its asymptotic value. Af is a scaling factor, and it is the asymptotic value which AfAt(t) approaches, as t gets large.
Here is the console output from the script, and the plot that it generates.
>> myModelTest
i=1, x=1.10e-10,0.900, rmsErr=0.0356
i=2, x=1.10e-10,0.960, rmsErr=0.0026
i=3, x=1.10e-10,1.000, rmsErr=0.0242
i=4, x=8.00e-11,0.960, rmsErr=0.0757
i=5, x=1.00e-09,0.960, rmsErr=0.3787
2 Comments
William Rose
on 22 Sep 2021
If you have a new set of experimental data which you want to fit, then save the numbers in a text file, one number per file, with no header, just like file AtAfExpt.txt. If the new file name is AtAfExpt2.txt, then modify fitData.m. Change
AtAfExp=load('AtAfExpt.txt');
to
AtAfExp=load('AtAfExpt2.txt');
I made one change in fitData.m and in fitDataPolyStart.m and myModel.m: I replaced "3001" with "length(AtAfEp)". That shoud assure a correct result, even if the new data file has a different number of numbers in it. That could have been the source of the incompatible size errors. Try the attached files, after you change the file name in fitData.m.
Categories
Find more on Large Files and Big Data in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
