Initial values in nlinfit or fitnlm

I am trying to run a non-linear multiple variable model in Matlab. The model has about 20 coefficients. I have been using 1s as my initial values in developing the model, and my model has an acceptable R2 value and good residual plots. However, I am not sure if the generated coefficients are sensitive to the initial values that were assigned by me.

This begs the question of whether one can check if the generated coefficient values are highly sensitive to the assigned initial values? Or one should do it manually, ie test the model with different set of initial values and compare the RMSE of the model?

Answer 1

[Star Strider]

A model with 20 parameters is likely going to be a challenge. If you have any doubts — and if you have the Global Optimization Toolbox — use the patternsearch (link) function to find the best parameter set. Another option is the genetic algorithm, the Global Optimization Toolbox ga (link) function.

I would also use the coefCI (link) function to determine if any of the confidence intervals for the coefficients (parameters) include zero, i.e. have opposite signs. If they do, they are not required for the model, since they are not statistically different from zero. This can help you ‘trim’ your model.

[wesleynotwise]

Hello Star Strider, nice to hear from you again!!! I've seen your reply hours ago, but I got distracted and couldn't reply to you quickly.

I have the Global Optimization Toolbox and I tried to run the code as below:

beta1 = ones (22,1)                 % initial value for 22 coefficients
x = patternsearch(modelfun1, beta1) % the patternsearch code
x = ga(modelfun1, 22)               % the ga code
mdl = fitnlm(tbl,modelfun1,beta1)   % fit into nonlinear

However, both ran into the same problem "Not enough input arguments", which is caused by "Failure in initial user-supplied fitness function evaluation. PATTERNSEARCH (or GA) cannot continue.", Just to let you know I have no problem when I run my model.

Can I check with you, if the above two codes work, does it mean I should use their coefficient values as my initial values in my regression model?

mdl = fitnlm(tbl,modelfun1,beta2)  
% beta2 is the results from GA or PATTERNSEARCH

Also, correct me if I'm wrong, I realise that the PATTERNSEARCH function requires the initial values, does that mean it may also giving me the same problem, i.e. the coefficients depend on what you assigned.

And, thank you for the suggestion for the coefCI function. I have it in my codes, as I wanted to use it to round the coefficient values to either 1 or 2 decimal points. I've never thought that it can also be used to trim my model! What a brilliant suggestion!!!!!

[Star Strider]

My pleasure.

Since you’re fitting your function to data, you have to introduce a cost function and minimise it with patternsearch or ga.

Example —

x = ...;                                    % Independent Variable
y = ...;                                    % Dependent Variable

RNCF = @(b) norm(y - modelfun1(b,x));       % Residual Norm Cost Function

where ‘b’ is your parameter vector.

Both functions will search as exhaustively as you let them to find the optimal parameter estimates.

You do not require the fitnlm function to estimate your parameters later, since the Global Optimization Toolbox functions that you decide to use will fit them about as well as can be expected. They are most likely to find the global optimum without your having to guess the initial values, so your fitnlm call will simply provide you with a model to use to present to coefCI.

You can round the coefficient estimates using the round function, to the number of decimal places you want. (I don’t remember when this option was introduced, so if your documentation for round doesn’t include it, I can post a one-line anonymous function that does the same thing.)

‘What a brilliant suggestion!!!!!’

Thank you!

[wesleynotwise]

Sorry, did I miss something in the code. I could not find where the patternsearch or ga function is in the example?

'You do not require the fitnlm function to estimate your parameters later' ---> Does that mean that fitnlm is not the function one should use when developing a model? I am little confused now, if patternsearch or ga function is better than fitnlm function in finding the parameter estimates, then what is the use of fitnlm?

I am using 2017a version, think I am okay with round function. Thanks:)

[Star Strider]

To clarify, I intend:

beta1 = ones (22,1)                % initial value for 22 coefficients
xp = patternsearch(RNCF, beta1)    % the patternsearch code
xg = ga(RNCF, 22)                  % the ga code
mdlp = fitnlm(tbl,modelfun1,xp)    % fit into nonlinear
mdlg = fitnlm(tbl,modelfun1,xg)    % fit into nonlinear

Use ‘mdlp’, ‘mdlg’, or both, with coefCI.

My pleasure.

[wesleynotwise]

Thanks for the clarification. I have some problems when I try to run the pattersearch or ga code, which is due to ' Undefined function or variable 'x' '. Not sure which part of my code went wrong?

beta1 = ones (22,1)                 % initial value for 22 coefficients
RNCF = @(b) norm(y - modelfun1(b,x))% Residual Norm Cost Function
xp = patternsearch(RNCF, beta1)     % the patternsearch code
xg = ga(RNCF, 22)                   % the ga code
mdlp = fitnlm(tbl,modelfun1,xp)     % fit into nonlinear
mdlg = fitnlm(tbl,modelfun1,xg)     % fit into nonlinear

Also, if the outcomes from patternsearch code and ga code are different, and so the parameter estimates of mdlp and mdlg, I assume I have to choose to use mdlp or mdlg based on its RMSE, as well as the R2 value?

[Star Strider]

You will have to adapt my code to your function and data.

Specifically note:

x = ...;                                    % Independent Variable
y = ...;                                    % Dependent Variable

That should clarify the ‘x’ and ‘y’ in my code example. I do not have your data or your function, so example code is the best I can do.

‘I assume I have to choose to use mdlp or mdlg based on its RMSE, as well as the R2 value?’

I would be surprised if they came up with the same parameter estimates. I would use the estimates with the lowest RMSE. The R² should track with the RMSE, so I doubt considering both would be helpful.

You can tweak the patternsearch and ga functions to search a bit more exhaustively. For patternsearch, that would be 'MaxIterations', and for ga, 'MaxGenerations'. Whether that would be beneficial depends on how quickly they converge. (The current default value would be 2200 for both.) With 22 parameters, they might need more iterations to arrive at the best parameter estimates.

[wesleynotwise]

Ah... no wonder. But, at the moment, I fit the model from a table but not a matrix. See the codes below. And the table that I built has more input than what I actually need for my model, as it is still pretty much in the development state. I assume I need to build a matrix for your code in order to incorporate in the existing one?

tbl = table(CN, CR, CON, FON, ANT, ART, AS, NAK, RSK, A_AN,...
              A_AR, ACa, SMPa, SAN, SAR, ANonly, CS, CNS,...
              EMO); % The table has more input than i actually need
modelfun1 = @(b,x)(((SMPa < 10).*b(1).*(x(:,13).^b(2))+...
                   (SMPa >= 10).* (x(:,13).^b(3))).*...
% the equation is very long, I only showed part of it   

X = [CN, CR, CON ...]; % I assume this is needed ?
y = EMO;               % And this?

beta1 = ones (22,1)                 % initial value for 22 coefficients
RNCF = @(b) norm(y - modelfun1(b,x))% Residual Norm Cost Function
xp = patternsearch(RNCF, beta1)     % the patternsearch code
xg = ga(RNCF, 22)                   % the ga code
mdlp = fitnlm(tbl,modelfun1,xp)     % fit into nonlinear
mdlg = fitnlm(tbl,modelfun1,xg)     % fit into nonlinear

Excuse the messiness in my codes. Need a good housekeeping.

[wesleynotwise]

Also, come to think of the patternsearch and ga functions, I think the latter is perhaps slightly better as it does not rely on the initial values, but just the number of parameter. Am I right?

[Star Strider]

x = [CN, CR, CON ...]; % I assume this is needed ?
y = EMO;               % And this?

Yes to both. Note that I corrected ‘X’ to ‘x’, since MATLAB is case-sensitive.

‘Am I right?’

Correct. The ga approach finds the best fit parameters from a very large initial parameter space, eventually coming up with the best overall set.

[wesleynotwise]

Yeah. Haven't run the codes, but I think it should work perfectly. Thank you very much! I really appreciate it!

Can I sidetrack a bit: All the three approaches, i.e. patternsearch, ga and fitnlm, have the same outcome which is the parameter estimates. We are now using the first two to find the 'suitable initial values' and use them in the fitnlm function. My questions are ->

Q1: Should these three functions do the same job, i.e. finding the optimum parameter estimates in a model? Or they have their own roles?

Q2: As you mentioned, ga approach starts from a very large initial parameter space, does it mean that the parameter estimates from ga is more 'correct' than fitnlm, given that the latter relies on the initial values?

[Star Strider]

My pleasure!

You can always sidetrack as you please. I will answer within the areas of my knowledge.

1. They have their own roles. The patternsearch and ga functions search the entire (or a very large part of the) parameter space for the best parameter estimates. The fitnlm function searches in the region near the initial estimates you’ve given it. The advantage of fitnlm is that it then allows you to calculate the statistics on the fit.
2. The parameters estimated by ga are more likely to be the most accurate, because it searches more widely. In a parameter space with a global minimum that is relatively straightforward to find, all parameter estimation routines will work optimally, and find essentially the same parameter estimates. The problem arises when there are several local minima that fitnlm, using a gradient-descent approach, could become ‘trapped’ in. Since ga does not use a gradient-descent approach, it is more likely to find the global minimum without getting trapped in local minima. When you then give those parameter estimates to fitnlm, it will converge quickly on the optimal parameter estimates, and give you the statistics on the fit.

[wesleynotwise]

Very very clear explanation. I am glad that I have this question in mind but now it is resolved beautifully. Thank you.

'within the areas of my knowledge' -> Oh well, I think your knowledge is as wide as the ocean.

[Star Strider]

Thank you. My pleasure.

I only wish I knew as much as I would like! The best I can do is to read and experiment.

[Priya Goel]

Hello Star Strider,

I went through your previous comments and followed your suggestions for my problem of curve fitting. I was able to implement patternsearch and ga but I faced some issues.

1. For patternsearch: The estimated parameters are sensitive to the initial guess. The RMSE value is same for multiple combinations which makes difficult to decide the best set of parameters.

1. For ga: It gives completely random set every time I run. The number of data points is only 1600. Can it be an issue for ga to not work properly? It also estimates negative values for some coefficients which is not logical as per the model. How can I define this constraint?

I will be really thankgul if you give some information on these.

[Star Strider]

1. Depending on the model, there may be several parameter combinations that provide essentially the same fit. This is likely a problem with the model.
2. It is usually necesary to run ga several times, ideally with a large initial population, to get the best fit. It is possible to automate this (running ga in a loop) and then saving the best individual from each run. See: How to save data from Genetic Algorithm in case MATLAB crashes? for the necessary code.

[Priya Goel]

Hello Star Strider,

Thank you for the quick reply. It is now clear that there may be more than one set of parameters which give same fit. But if this issue arises because of the model, it will continue to exist even if ga is used.

[Star Strider]

My pleasure.

Very true! That is usually because the model is not well formulated, and does not correctly represent the process that created the data.

[Priya Goel]

In continuation to your statement, can it also happen when data is not sufficient. To elaborate, I have two independent variables in the model. Now, the data has two type of entries:

a) First independent variable changes and second is constant

b) Second independent variable changes and first is constant

Number of entries for case a) is sufficient but it is very less for case b). It may lead to poor fit for the coefficient of second variable.

[Star Strider]

Pardon me, but that does not make any sense!

If you have more data than parameters (preferably at least twice as many data as parameters), the data are not excessively corrupted by noise, and the model acccurately represents the process that created the data, there should be a decent fit.

[Priya Goel]

i am sorry if this seems to be unreasonable. But I want to be more clear in this. Let us suppose we are trying to fit data using an equation:

Z = A*(x^b)*(y^c)

Now, there are 3 coefficients which need to fit. In order to get a decent fit, I should have data that contains variation of x and y. However, if y is constant for most of the entries, I may not be able to obtain c accurately.

[Star Strider]

This should be a new Question.

It does not directly relate to the current thread. I am not going to respond to it further here.

[Priya Goel]

Yes. I am also of the opinion that it is deviation from the original topic.

Anyways, Thank you for your valuable time and inputs. These are very helpful for beginners (like me). You clarify doubts which otherwise remain unanswered for months.