set parameters of nlinfit function

Hi

I’m beginner in MATLAB. I have some experimental points and I used nlinfit function to solve my problem and find the best LS curve fitting result. When I use the nlinfit, it pass MSE= 1.5873e+07 and the result is not very good. I attached my data, please check output of nlinfit function to see the output.

Note: I want to find the best curve fitting based on spherical model, soI use: fun = @(p,X) p(1)*(1.5*X/p(2)-0.5*(X/p(2)).^3);

How I should set function’s arguments to reach the best results? Please help me to solve my problem.

Best regards

Answer 1

[Star Strider]

One problem might be that you did not completely vectorise your objective function code (notably the divisions and at least one multiplication. I did here, so that could help:

fun = @(p,X) p(1).*(1.5*X./p(2)-0.5*(X./p(2)).^3);

See if running it makes a difference. If it solves your problem, let me know.

I’d actually like to know a bit more about what you did. It would help to know your starting parameter estimates.

Past midnight here, so I’ll run your code with your data later in the morning and see what results I get.

[Mani Ahmadian]

Dear Star Strider;

thanks a lot for your note about my mistake. I correct it, but nothing improved. I attached my function, too. Please run it with my own data and tell your idea.

Sweet dreams... Mani

[Star Strider]

This code works, as I assume yours does:

fidi = fopen('nlinfit_Data.txt');
D = textscan(fidi, '%f %f', 'HeaderLines',1, 'Delimiter',' ' );
xdata = D{:,1};
ydata = D{:,2};

fun = @(p,xdata) p(1)*(1.5*xdata./p(2)-0.5*(xdata./p(2)).^3);
pguess = [1; 2];
b = nlinfit(xdata,ydata,fun,pguess);

yhat = fun(b,xdata);

figure(1)
plot(xdata, ydata, 'xb')
hold on
plot(xdata, yhat, '-r')
hold off
grid

You may not have chosen the correct model, but since I have no idea what your data are or what you want to do with them, I can only suggest you explore other models.

[Mani Ahmadian]

Thanks a lot for your guid.

All of points are results of variogram calculation. I have to choose one of theoretical model as Spherical, Exponential or Gaussian.

Each above models have specific function as: model's function=f(h,a) %mathematical expression

Spherical model:

Exponential model:

Gaussian model:

I changed your code as attached file. After run code, I find 2 Points: 1- I have some warnings in Gaussian and exponential models, but I don't know how I can remove them. 2- Two last models are sameness(I guess results are not true because of above warnings).

I think spherical model is the best one based on diagrams in plots. Is it true?

Another item is to change parameters of nlinfit function to move up start point of estimated diagram to close it to y=0.5*10^4. Is it possible?

[Star Strider]

My pleasure!

‘I think spherical model is the best one based on diagrams in plots. Is it true?’

Don’t depend on how the plots look. This is a statistical exploration.

Note that nonlinear parameter estimation routines are sensitive to the initial parameter estimates.

Use:

pguess = [2E+4; 1E+3];

for your initial estimates for every model, and you may conclude that other models may offer closer fits to your data than your earlier initial estimates provided. (The exponential fit starts looking better with those initial estimates.)

I also suggest you return to your original syntax in your call to nlinfit:

[p,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(xdata,ydata,fun,pguess)

because you need those outputs to input to nlparci and nlpredci and other calculations to help you determine the best fit. Remember to name them slightly differently in each call to nlinfit so you don’t overwrite earlier results in your subsequent calls to nlinfit.

As to the best fit, that is entirely up to you to determine. (The norm of the residuals is a good place to start.) They are your data and regressions and statistics.

[Mani Ahmadian]

Great... I changed the initial values by yours and I see the wonderful effect. Thank you very much.Please guide me how I guess the initial values, specially about spherical model(red line) because it doesn't change and I have to improve it and then choose the best one.

As another point, please introduce me a beneficial text book to improve my knowledge about "Norm of residual", I have to learn more and more...

Best regards, Mani

[Star Strider]

My pleasure!

‘Guessing the initial values’ involves knowing something about your data and the function you are fitting. In your application here, the x- and y-data are on the order of 1E+4, so it makes sense to initialise the parameters to a similar magnitude, considering that your parameters are scaled by the x-data in all your functions. (For instance, if you multiplied your data by your parameters instead of dividing your data by the parameters, the appropriate initial parameter estimates magnitudes would be on the order of 1E-4 instead.) The important considerations in choosing initial parameter estimates for a particular problem are understanding the maths of your regression objective functions and the data you are using in your nonlinear parameter estimations.

There are a number of good textbooks on linear and nonlinear parameter estimation and regression that use MATLAB. See: MATLAB and Simulink Based Books. There are a number of online resources available as well. It is probably best that you search for those yourself, since you have a much better understanding of your background and applications than I do, and while having an encyclopaedic knowledge of these techniques is desirable, it is not usually practical.

Another source I use to guide me to good references are university courses on the subjects I’m interested in. If the course looks as though it covers the topics that interest me, I search for the syllabus. The professors usually mention the textbooks for the course in the syllabus for the course.

[Mani Ahmadian]

Dear Star Strider

Thank you very much for your useful guidance, I'll try to search and study more to improve my statistical knowledge. I'm so pleased to meet you here.

About my problem, by use of your nice idea I could solve pitfalls about Exponential and Gaussian models, but spherical model does not change. Would you please help me to complete it,too?

Have a nice time, Mani.

[Star Strider]

My pleasure!

If you haven’t already done so, I suggest you start with a good introductory Statistics course. Then depending on what you want to do, take specialty courses, for instance in linear and nonlinear regression if you want to continue with what you are doing now.

The Spherical model seems to fit its parameters well even if given an initial estimate that is far from ideal. (You would have to calculate its Jacobian to see the reason.) I used various initial estimates for it, and it always converged on the same values. Assume it is already optimised as much as it can be.

I would use the norm of the residuals or analogously ‘MSE’ (mean squared error) as a way to determine the best model, so the ‘best’ model would have the least MSE. (If you want, you can estimate the statistical significance of the difference between two models with the Likelihood Ratio Test. For your purposes, it is safe to assume the the MSE is the likelihood, since all your models have the same number of parameters. This is the only way I know to compare two nonlinear models with different objective functions.)

Did we solve your problem or do you still have questions?

[Mani Ahmadian]

Of course... you helped me to solve the most parts of my problem. But condition of spherical model doesn't satisfied me. I think it's diagram is out of tune and its mse is a huge value. I guess I have a big calculation or mathematical mistake about it, but I'm not sure where probably it happens.

I downloaded some textbooks about nonlinear curve fitting and I'll study them as soon as possible to have better understand.

Your view points and guidance was very helpful for me. Thank you very much, my friend.

Have a good time. Mani

[Star Strider]

My pleasure!

Compare the MSE values for the various models. If any of the models adequately explain the process that produced your data, those should take precedence. If they all do, then the one with the lowest MSE is likely the correct one. If with the Likelihood Ratio Test they are not significantly different from the one with the lowest MSE, they are all valid. Developing and estimating other models that could explain your data are appropriate and acceptable. If you can find a more appropriate model, fit it as well and compare it to the others. That is the fun — and frustration — of nonlinear parameter estimation and mathematical modeling.

You have fun as well! I wish you well in your research. We’re here if you need us for other MATLAB projects.

Answer 2

[Mani Ahmadian]

Hi again

I'm trying to solve the problem with spherical model fit function, this model has a mathematical formula as below:

I think the un-fit problem is due to two criteria function of above model. I don't know how I should define this type of function to nlinfit.

Would you please help me to solve the problem?

Thanks

Mani

[Star Strider]

I would estimate it only for h <= a0, since that estimates c0 as well.

Defining the parameters as:

% b(1) = c0,  b(2) = a0
gam2 = @(b,h) b(1).*(1.5*h./b(2) - 0.5*(h./b(2)).^3);

with the function ‘gam2’ being the function you would use as your objective function for nlinfit.

[Mani Ahmadian]

Dear Star

By use of your great comments above I learned how to use nlinfit function. But now, I need to define the model function as a two critical function.

About your idea in last comment, I think it is not sufficient, because nlinfit fun. never knows about the second criteria, which is a horizontal line after h>a0. So I think I should change my function definition to nlinfit to have the best fit.

As I told in my previous comments, about Gaussian and exponential, nlinfit function has the best fit. I guess nice result of this function is due to one criteria function of them. I used same data on all models and test many ways to improve nlinfit result about spherical model, but nothing changed!

Now I'm seeking to find a way or any other mathematical ways to find best fit for spherical model.

Please help me to solve the problem.

Thanks

Mani

[Star Strider]

I continue to suggest that you fit only the values for h <= a0 because of the discontinuity at h = a0. The fitting functions must be continuously differentiable, since nlinfit uses a Levenberg-Marquardt gradient descent algorithm, and they are not at such discontinuities. You can estimate c0=b(1) as well as b(2)=a0 on the continuous portion of the curve.

[Mani Ahmadian]

Dear Star

I'm so sorry for my late answer, I had a trip last week. About your answer, I think it's not a complete answer because if we have a function, y, as bellow:

As we discussed above, I write a simple script to use nlinfit function. I use nlinfit function in two steps, in first step, I use y1=x^2+x+4 and for second step I use y2=abs(x)+4

As you see, the result of nlinfit function is not perfect while we are using just one of y1 or y2. So I think to reach the best result, we should find a way to input two function of y1 and y2 simultaneously.

Please check my script and tell your idea.

Thanks

Mani

[Star Strider]

You cannot fit a continuous and continuously-differentiable function across a discontinuity such as you describe. You have to fit each region of the discontinuous function separately, possibly with a different function describing each continuous region.

That is just how the maths work. It is beyond my power to change them!

[Mani Ahmadian]

Thanks for your fast answer

Please consider the spherical model formula in above answers. As we discussed, finding x=a is the most important pitfall.

I have an idea. Please imagine we have 20 data points, because for x<=a the spherical function doesn't have an constant value, I want to use a loop and nlinfit with 3,4,5,6,...,20 data points. In each stage I collect the norm of the residuals,...

Finally, I want to select the best result of norm of the residuals. In this way, I will find value of a and for x>a we can use constant value for spherical model.

What's your idea about my approach?

[Star Strider]

It is best to fit all of the data at once.

You set a threshold of h<=a0 but haven’t defined what ‘a0’ is. Fitting it then is simply a matter of using only those (h,gamma) data.

What is ‘a0’?

[Mani Ahmadian]

a0 is x value of a point and for points with greater that a0, y is constant and is equal to C0.

[Mani Ahmadian]

I should find a way to read a0 from data points and result of nlinfit function. My reason to use loop is to find x=a0 so that for x<=a0, the spherical model has an non-linear formula.

[Star Strider]

Truncating the x and y values so that x<=a0 for successive values of ‘a0’, both parameter values converge to essentially stable values, although those values are a function of ‘a0’.

The following code (for the spherical model) loops, calculating ‘a0’ at each step and truncating the (x,y) arrays to accommodate x<=a0 in the subsequent step, estimating the parameters with the new data set, and repeating for 25 iterations:

fidi = fopen('nlinfit_Data.txt');
D = textscan(fidi, '%f %f', 'HeaderLines',1);
x = D{1};
y = D{2};

% b(1) = c0,  b(2) = a0
gam2 = @(b,h) b(1).*(1.5*h./b(2) - 0.5*(h./b(2)).^3);

B(:,1) = [1E+4; x(length(x))];
Lxy(1) = length(x);
for k1 = 2:25
    B(:,k1) = nlinfit(x(x<=(B(2,k1-1))),y(x<=(B(2,k1-1))),gam2, B(:,k1-1));
    Lxy(k1) = length(x(x<=(B(2,k1-1))));
end

figure(1)
semilogy(1:25, B(1,:),  1:25, B(2,:))
legend('c_0', 'a_0')
xlabel('Iteration #')
ylabel('Parameter Value')

The ‘B’ matrix contains the parameter values for each iteration, and the ‘Lxy’ vector (length of the data vectors at each step) demonstrates that the data are indeed being truncated. The preceding parameter values are used as the initial estimates for the subsequent estimation. The results are yours to interpret.

To plot the estimated y-values for a particular index of parameter values, you can use this code:

Lidx = 9;                           % Choose Parameter Set By Index
xv = linspace(min(x),x(Lidx));      % Create x-Vector

yfit = gam2(B(:,Lidx), xv);         % Calculate Estimated Fit

figure(2)
plot(x(1:Lidx), y(1:Lidx), '+')
hold on
plot(xv, yfit, '-r')
hold off
grid
xlabel('\ith\rm')
ylabel('\gamma_{2}(\ith\rm)')

Have fun with it!

[Mani Ahmadian]

Great

Thanks for your nice job. I'll check it as soon as possible.

Have a nice time

Mani

[Star Strider]

My pleasure!

You too!