Curve fitting with custom function

Hello community,

I am trying to do a curve fitting on some experimental data with a custom function. I am more specifically trying to bring closer my function to the data. The function is the following one:

Every parameters of the function, except the shear rate () , can vary in order to fit best the data.

Here is a graph with the data and the function plotted with initial values:

The initial values are the following ones:

muInf=0.0035 %[Pa.s]

mu0=0.108; %[Pa.s]

lambda=8.2;

n=0.3;

a=0.64;

viscCar=muInf+(mu0-muInf)*(1+(lambda.*shearRate).^a).^((n-1)/a);

loglog(shearRate,viscCar,'k');

I already tried to use the curve fitting toolbox but i wasn't able to keep the look of the function.

Does anyone know how i can do that ?

Thank you for the help !

Answer 1

[Star Strider]

It would help to have the data.

I would use the fitnlm function for this —

shearRate = logspace(-3, 4, 100);                                   % Create Data
muv = 0.5-tanh(shearRate)*0.01 + randn(size(shearRate))*0.01;       % Create Data

shearRate = shearRate(:);
muv = muv(:);

% muInf=0.0035; %[Pa.s]
% mu0=0.108; %[Pa.s]
% lambda=8.2;
% n=0.3;
% a=0.64;

B0 = [0.0035; 0.108; 8.2; 0.3; 0.64];

viscCar = @(muInf,mu0,lambda,a,n,shearRate) muInf+(mu0-muInf)./(1+(lambda.*shearRate).^a).^((n-1)/a);
viscCarfcn = @(b,shearRate) viscCar(b(1),b(2),b(3),b(4),b(5),shearRate);

mumdl = fitnlm(shearRate,muv, viscCarfcn, B0)

Warning: Rank deficient, rank = 4, tol = 1.142224e+00.

Warning: Rank deficient, rank = 4, tol = 1.739422e+00.

Warning: Rank deficient, rank = 4, tol = 2.345437e+00.

Warning: Rank deficient, rank = 4, tol = 3.085596e+00.

Warning: Rank deficient, rank = 4, tol = 3.036782e+00.

Warning: Rank deficient, rank = 4, tol = 4.094796e+00.

Warning: Rank deficient, rank = 4, tol = 5.447783e+00.

Warning: Rank deficient, rank = 3, tol = 5.383975e+00.

Warning: Rank deficient, rank = 1, tol = 6.115407e+01.

Warning: Some columns of the Jacobian are effectively zero at the solution, indicating that the model is insensitive to some of its parameters. That may be because those parameters are not present in the model, or otherwise do not affect the predicted values. It may also be due to numerical underflow in the model function, which can sometimes be avoided by choosing better initial parameter values, or by rescaling or recentering. Parameter estimates may be unreliable.

mumdl = 
Nonlinear regression model:
    y ~ viscCar(b1,b2,b3,b4,b5,shearRate)

Estimated Coefficients:
          Estimate         SE          tStat        pValue 
          _________    __________    __________    ________

    b1      0.23921    3.7712e-33    6.3429e+31           0
    b2      0.45235    2.3787e-32    1.9017e+31           0
    b3    1.391e-17    5.3563e-18        2.5969    0.010839
    b4     0.036323    5.3121e-32    6.8378e+29           0
    b5      0.97281    3.2215e-32    3.0197e+31           0


Number of observations: 100, Error degrees of freedom: 99
Root Mean Squared Error: 0.014
R-Squared: -0.682,  Adjusted R-Squared -0.682
F-statistic vs. zero model: 1.25e+05, p-value = 2.2e-155

Beta = mumdl.Coefficients.Estimate

Beta = 5×1

0.2392 0.4523 0.0000 0.0363 0.9728

figure
loglog(shearRate,viscCarfcn(Beta,shearRate),'k');
hold on
plot(shearRate, muv, '.b')
hold off
grid

This works, and would actually make sense with the actual data.

.

[Da125]

Thank you Star Strider for your answer !

Yes i attach the .csv file with the data.

I tried to implement your suggestion but i have issues with that. I obtained some negative coefficient that don't make sense. Here is what i obtained:

With the coefficients:

[Star Strider]

There were several errors in ‘viscCar’ that I found and coirrected. The most important was replacing the vectorised multiplication of the two terms by a vectorised division, and the others were adding appropriate parentheses to the denominator and replacing ‘(n-1)’ with ‘(1-n)’ in the expoenet. With those corrections, it works. (I added the fminsearch call to estimate ‘B0’ for fitnlm because fminsearch uses a derivative-free search, and would likely find a suitable set of parameters relatively quicklly that it can then provide to fitnlm.)

The principal problem, however, was mis-coding the ‘visccar’ objective function.

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/803909/bloodExpData.csv', 'VariableNamingRule','preserve');

shearRate = T1.Var1;
muv = T1.Var2;
% shearRate = logspace(-3, 4, 100);                                   % Create Data
% muv = 0.5-tanh(shearRate)*0.01 + randn(size(shearRate))*0.01;       % Create Data

shearRate = shearRate(:);
muv = muv(:);

% muInf=0.0035; %[Pa.s]
% mu0=0.108; %[Pa.s]
% lambda=8.2;
% n=0.3;
% a=0.64;

% B0 = [0.0035; 0.108; 8.2; 0.3; 0.64];
B0 = rand(1,5);

viscCar = @(muInf,mu0,lambda,a,n,shearRate) muInf+(mu0-muInf)./((1+(lambda.*shearRate).^a).^((1-n)/a));
viscCarfcn = @(b,shearRate) viscCar(b(1),b(2),b(3),b(4),b(5),shearRate);

[B0,fv] = fminsearch(@(b) norm(viscCarfcn(b,shearRate) - muv), B0)

B0 = 1×5

0.0039 0.6823 0.1168 0.1792 -0.3050

fv = 0.0019

mumdl = fitnlm(shearRate,muv, viscCarfcn, B0)

Warning: Iteration limit exceeded. Returning results from final iteration.

Warning: The Jacobian at the solution is ill-conditioned, and some model parameters may not be estimated well (they are not identifiable). Use caution in making predictions.

mumdl = 
Nonlinear regression model:
    y ~ viscCar(b1,b2,b3,b4,b5,shearRate)

Estimated Coefficients:
          Estimate         SE          tStat         pValue   
          _________    __________    __________    ___________

    b1    0.0035412    0.00017069        20.746      2.846e-17
    b2       462.39    8.1834e-07    5.6503e+08    7.4334e-203
    b3        789.3    7.8721e-07    1.0027e+09     4.407e-209
    b4     0.083341     0.0029383        28.363     1.5629e-20
    b5      0.15071      0.020071        7.5092     7.3063e-08


Number of observations: 28, Error degrees of freedom: 25
Root Mean Squared Error: 0.000347
R-Squared: 0.997,  Adjusted R-Squared 0.996
F-statistic vs. zero model: 9.08e+03, p-value = 4.37e-38

Beta = mumdl.Coefficients.Estimate

Beta = 5×1

0.0035 462.3854 789.3003 0.0833 0.1507

figure
loglog(shearRate,viscCarfcn(Beta,shearRate),'k');
hold on
plot(shearRate, muv, '.b')
hold off
grid

Success!

.

[Da125]

It works ! thank you !

[Star Strider]

As always, my pleasure!

.

Answer 2

[Alex Sha]

It is hard to get stable and unique result for Da125's problem, especially for parameters of "lambda" and "n". refer to the result below:

Root of Mean Square Error (RMSE): 0.00032178294087402
Sum of Squared Residual: 2.89923930905092E-6
Correlation Coef. (R): 0.998376698597668
R-Square: 0.996756032302778

Parameter	Best Estimate
----------	-------------
muinf	-0.0902678665303079
mu0	0.00320052355322768
lambda	207505339.640542
a	-0.504386079091155
n	-1252.71593136491

[Star Strider]

Negative parameters are apparently not permittted. This was a problem noted in a previous Comment.

.