How to fit curves with two variables to equation?

Hello Matlab-Community,

I am struggling quite a while with the problem of getting a formula out of 3 graphs I obtained from simulation.

Below you see the quantity k over normalized concentration for three different velocities.

The velocities vary in the range of 0.001, 0.01 and 0.1.

Using the fit function in matlab, I get a good approximation for the respective curves regarding k in dependency of c/c0.

f1=fit(c/c0,k,fourier5)

Now, I want to express k as a function of c/c0 and v. k = f(c/c0,v)

I checked several threads but could not manage to obtain an equation.

My problem is also that the three curves are built up of different sized matrices (165x3 double, 80x3 double and 83x3 double) with different interval size.

Do you have any hints?

Best Rgards,

Amadeus

Answer 1

[William Rose]

@Amadeus Wolf, please post the data (the three matrices) which you are trying to fit: k and c/c0 and v.

[edit: Remove reference to fmincon(). Add mention of custon function as a fitType. Fix typos.]

On the plot, I recommend distinguishing between the data and the fit.

If you want to use the fit function, then you are fitting a surrface: k is a function of two variables, c/c0 and k. Therefore you cannot use the fourier5 model, which is for functions of one variable only. You can use lowess or loess or polymn models, where m=1 to 5 and n=1 to 5:

f1=fit([c/c0,v],k,fitType)

where fitType is replaced with loess, lowess, poly11, poly 21, poly12, ..., or poly55.

The diffrent array sizes are not a problem. You can combine your three arrays into a single array with size 245x3, whose columns correspond to c/c0, k, and v.

If none of these fitTypes give a good fit to your three curves, you may have to make your own custom fitType function, and still use fit() to fit it. See here.

[William Rose]

@Amadeus Wolf, Inspect the fitted coefficients for each of the v-values which you have fitted with fourier5. There should be 11 fitted coefficients with Fourier5: a0, a1, ..., a5 and b1, ..., b5. See if you notice a pattern with respect to the velocities. THis can help you write your own fitting function.

Your data appears to be symmetric, or very close to symmetric, about 0.5, on the interval c/c0=0 to 1. If that is the case, then the sine terms in the Fourier fit should all be zero. Consistent with this, I predict the values you have already found for b1, ...b5 will be small. Therefore, if you make your own function to fit the surface, I recommend not including sine terms. Maybe something as simple as the following will work well:

or

where x stands for c/c0.

[William Rose]

@Amadeus Wolf,

I made a text file with 36 data points (36 triples of v,c/c0,k) that resemble the points on your blue, red, and black curves.

I fitted the function below to the data. This function has only 4 adjustable parameters. You fitted one curve at a time with fourier5, which has 11 adjustable parameters. A model with fewer adjustable parameters is preferable, other things being equal.

where a,b,c,d are the parameters that get adjusted in order to get a good fit, and , which is called cc0 in the code, is the horizontal-axis-value on your original plot and on the plot below. I chose the functional form above because I thought it might fit pretty well, with a small number of adjustable paramaters.

The code to do the fit is below.

a=load('k_c_v.txt');
v=a(:,1); cc0=a(:,2); k=a(:,3);
Effect = @(a,b,c,d,x,y) c+d*x+b./sin(pi*a*y);   %custom fitting function; x,y correspond to v,cc0
[f,gof] = fit( [v,cc0], k, Effect, ...         %fit: custom
    'StartPoint', [1, 0.002, 0.01, 0.02], ...
    'Lower', [0,0,-0.1,-1],'Upper',[3,3,0.1,1]);
kpred=feval(f,[v,cc0]);                         %predicted values with fit
fprintf('Custom Fit: rmse=%.5f.\n',gof.rmse)    %display root mean square error

Custom Fit: rmse=0.00126.

The resulting fit is shown in the plot below. The X's are the data points we are trying to fit, and the lines are the fit. A single model with four parameters is used to fit all the data.

[Amadeus Wolf]

Hi @William Rose,

thank you very much for the detailed response.

I also tried to reduce the number of parameters of my fit. Using poly22 or higher in the surface fit did not lead to a good result.

Regarding the 2D fit, a Holliday-type function lead to an appropriate result so far:

Unfortunately, this function type does not work for surface fit:

"Custom equations must produce an output vector, matrix, or array that is the same size and shape as the input data. This custom equation fails to meet that requirement: "

Do you know, why this happens?

I attached some manipulated fake data of mine. Hope that is also fine.

I will check out, if your function will fit the data. Is it vaild for the intervall c/c0 [0,1] ?

Again, thanks al lot!

Best Regards,

Amadeus

[Amadeus Wolf]

@William Rose Addition:

I want to extrapolate the curves for the intervall c/c0 [0,1]. For the missing range, the plan is to use the values from the fourier5 fit. They fit the best.

Using this "extended" curves I want to obtain the final formula k = f(c/c0,v) in the whole intervall c/c0 [0,1] and v [v1,v3] to use it somewhere else.

[Amadeus Wolf]

I tried to use your approach on my data but I it did not worked out.

I think it is because of the Start Point, upper and lower limit. On what basis did you define them?

Attached you find the final and extrapolated curves in a file.

Best Regards.

[William Rose]

@Amadeus Wolf,

I determined the start point for the fit by looking at the plotted data in your original posting, and experimenting in Excel. The data in your plot were for v=1e-3, 1e-2, 1e-1. The data in file k_c_v_AW.txt are for v=5e-5,5e-4, and 5e-3. Also, the data in the orignal plot did not include c/c0=0 or c/c0=1. The function I have proposed goes to +Inf when c/c0=0. This was not a problem with the orignal plot data but it is a prblem with the mroe recent data you have posted. These dufferences between the newer data and the orignal plot are why the initial guess does not work for the data in the text file.

Please explain the relationship between the four files you have attached in the comments above. Is file k_c_v_AW.txt the combination of files cc0_k_v1_fake.txt, cc0_k_v2_fake.txt, cc0_k_v3_fake.txt? If not, then how does k_c_v_AW.txt relate to the *fake.txt files? Which files contain the simulated data which you want to fit, and which contain data that is the output from a fitting attempt?

You asked:

Unfortunately, this function type does not work for surface fit:

"Custom equations must produce an output vector, matrix, or array that is the same size and shape as the input data. This custom equation fails to meet that requirement: "

Do you know, why this happens?

I don't know why it happens, because you have not shown the code you used to implement the function. I suspect that it happens because the dot operator is not used in the necessary positions. If x and v are n-by-1 vectors, and if a1, a2, a3, and b1 are scalars,

x=(0:.05:1)'; v=.005*ones(21,1); a1=1; a2=1; a3=1; b1=1;

then

kpred = a1./(1+a1*x+a2*x.^2) + b1*v; disp(size(kpred))

    21     1

will also be n-by-1, and you should not get the error above.

[William Rose]

@Amadeus Wolf,

I have inspected the files cc0_k_v1_fake.txt, etc. It appears that the data in column 1 are not c/c0, but rather 100*c/c0. Is that correct?

I think the middle column has the k values. There are only 4 digits after the ".", and this means there are only 2 significant digits for many of the k values. I recommend including more digits in the column for k.

[William Rose]

@Amadeus Wolf,

The simulated data values do not include c/c0=0 or c/c0=1. The k values in tthe simulated data go up sharply as c/c0 approaches 0 and 1. Is there a theoretical reason why k would or would not be infinite or finite at c/c0=0 and 1? I know that the Fourier extrapolation predicts finite values, but a) extrapolation is often unreliable and b) just because the Fourier expansion predicts finite values does not mean the actual modeled phenomena would be finite at the limits.

I am asking because we want to chose a fitting function whose limiting behavior is consistent with the expected limiting behavior.

[William Rose]

@Amadeus Wolf,

I notice that the values for k in k_c_v_AW.txt are almost exactly 10 times smaller than the corresponding values for k, at the same c and v, that appear in cc0_k_v*_fake.csv, (*=1,2,3). I assume the k values in one of the files are wrong - which?

Also, as I mentioned before, the cc0 values in cc0_k_v*_fake.csv are multipled by 100 and rounded to the nearest imteger. As a result, c/c0 values of 91, 92, 93, 95, 96, 97 are all repeated in the file cc0_k_v1_fake.csv, and other c/c0 values in other cc0_k_v*_fake.csv files are also repeated. I am sure the actual c/c0 values were not repeated, but rounding makes them look the same. This will mess up the fit.

[William Rose]

@Amadeus Wolf,

I took your file k_c_v_AW.txt and I removed the values that wer obtaibed by Fourier5 extrapolation, because I don;t want to fit data that came from another fit. I only want to fit data from a simulation or from an experiment. The modified file is k_c_v_AWne.txt. (ne for no extraploated data.)

I fitted it with the same model as before, but with a different starting point. The diferent starting point was necessary because the k values and the v values are different by at least an order of magnitude from the values in the initial posted plot.

This resulted in a reasonable custom function fit with four adjustable parameters. The console output is below.

>> kcvFit

Fit 1 (poly14 ): rmse=4.30e-05.

Fit 2 (custom1): rmse=1.54e-05.

Fit 2 parameters: a=1.00147, b=4.09e-05, c=1.71e-04, d=8.25e-03

The plots produced are below.

The data file and script are attached.

[Amadeus Wolf]

Hi @William Rose,

thanks for your effort and the code you provided, I really appriciate it.

Regarding the function not working for the surface fit, you answered:

I suspect that it happens because the dot operator is not used in the necessary positions.

That's right. The code is

g2 = fittype(@(a0, a1, a2, b1, x, y) a0./(1+a1*x+a2*x.^2)+b1*y, ... 
     'independent', {'x', 'y'},'dependent', 'z' );
f1=fit([cc0,v],k,g2)

Forget the files cc0_k_v1_fake.txt, etc., only the k_c_v_AW.txt matters ;-)

Regarding the code you provided, it does the job pretty well. Thanks.