How to fit experimental data to a model defined by m. file

Hello everyone,

I am trying to fit experimental data to a non-linear model with two independent variables (z-position and t-time) and one dependent variable (C-concentration). So, my experimental data is defined as C=f(z,t), where each pair of value (z, t) leads to a unique value of C(z,t). I have 1000 z-values and 250 t-values. I have tried to tackle this task by using nlinfit (beta=nlinfit(X,Y,modelfun,beta0)), but it only works when X and modelfun are column vectors. Then, I combined the two independent variables (z, and t) into a matrix, and after that in a column vector as follow:

[t,z] = meshgrid(t,z);
X = [z(:); t(:)];

Now, the next step is to define modelfun. My model is quite complex to define it as a function handle, so I decided to write it in a separate file (fitDiff_nlinfit.m). The model function is:

where t_0 is an empirical parameter (to be regressed), t_end is the last time instance, and f is defined by:

H, b are constant for a given problem, and I will consider only the first 100 term in n.

Here, I would like to regress three parameters . So far, I have the following code:

D = 1e-09;
t_0 = 10;
DeltaC_0 = 1e-03;
varNames = {'D';'t_0';'DeltaC_0'}; 
fitPars = [D t_0 DeltaC_0];

[beta,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(X, C(:), ...
         @(beta)fitDiff_nlinfit(beta,varNames,t_0,t,z,H,b,D,DeltaC_0), fitPars);

But I received the error:

Error using nlinfit (line 213)
Error evaluating model function '@(beta)fitDiff_nlinfit(beta,varNames,t_0,t,z,H,b,D,DeltaC_0)'.

Error in calcDiff_nlinfit (line 126)
[beta,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(X, C(:),@(beta)fitDiff_nlinfit(beta,varNames,t_0,t,z,H,b,D,DeltaC_0), fitPars);

Caused by:
    Error using calcDiff_nlinfit>@(beta)fitDiff_nlinfit(beta,varNames,t_0,t,z,H,b,D,DeltaC_0)
    Too many input arguments.
    

what am I doing wrong? I attached here the scripts with the model function.

I appreciate in advance any help you can give me.

[dpb]

The function definition for the target for nlinfit can have only two arguments...the coefficient vector of coefficients being estimated and the array of X predictor variables. You've got a whole slew of arguments in your function, not just the two which triggers the error.

So, you need a footprint like

function yhat=myFunctionModel(beta,x)

that takes the coefficient vector beta and the independent variables array x and returns the estimate response as a column vector of response(s). In your case if I understand you only have the one response value but two predictor variables.

You must get other parameters into the functional by way of one of the techniques shown at PassingExtraParameters Unfortunately, the Stat TB documentation is awfully weak on its examples to show this issue; the above doc link is in the Optimization TB instead.

I've got to run at the moment so can't see if I can figure out your code sufficiently to recast it right now; if it's possible to do so, I'd suggest attaching a .mat file of your data too, somebody else may come along and be intrigued and attack it before I have a chance to get back...

But, that's the error you have now--your function footprint doesn't match the required form.

[Ronny Rives]

Thank you dpb, you're rigth about the lack of documentation on this problem. Now I will try to rearrange my model function. In any case, here I attach a .mat file with experimental data as a test in case someone is interested in attacking this example.

Answer 1

[Ronny Rives]

Hello,

I have made some changes in my code as suggested by dpb, and it`s working.

D = 1e-09;
t_0 = 10;
DeltaC_0 = 1e-03;
varNames = {'D';'t_0';'DeltaC_0'}; 
fitPars = [D t_0 DeltaC_0];

[t,z] = meshgrid(t,z);
X = [z(:); t(:)];

f = @(beta,X)fitDiff_nlinfit(beta,varNames,X,H,b,D,t_0,DeltaC_0);

[beta,R,J,CovB,MSE,ErrorModelInfo] = nlinfit(X,C(:),f,fitPars);

However, now I got the following warning:

Warning: Rank deficient, rank = 2, tol =  4.731182e-03. 
> In nlinfit>LMfit (line 587)
  In nlinfit (line 284)
  In calcDiffusionLong_nlinfit (line 132)
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.

One of the three parameters (in this case, t_0) remain invariable after the regression. The values of the other two regressed parameters (D and DeltaC_0) are physically consistent. However, when I changed the initial guess of t_0, the values of (D and DeltaC_0) also change, indicating that the model is sensitive to t_0.

Then, I tried to find better initial guesses as recommended. To do this, I used the fminsearch algorithm. The estimated parameters using fminsearch were physically consistent. So, I could keep these parameters values but I am interested in determining the standard deviation of the parameters. I know that it is possible to obtain the standar deviation of the estimated parameters with nlinfit (from CovB), but I don`t know if it is possible with fminsearch?

If it is not possible to obtain the standard deviation with the fminsearch algorithm, what else could I do to regress t_0?

[dpb]

I had returned to see if could help recast the function but see you cleared that hurdle, apparently. But, your attach of the .mat file didn't happen so can't see the data to judge what might be going on numerically.

The first suggestion w/o anything else would be to try one of the suggestions of the message -- use the estimates from fminsearch as your initial guess first, second recast the model to standardize the independent variables. How well either might work is impossible to guess without the data, of course.

[Ronny Rives]

Hi dpb,

Here you can find some experimental data. I tried to upload the file before, but it seems that it takes some time, and I closed the page before the upload was finished. The parameters from fminsearch are:

D = 7.573e-10;

t_0 = 45.84;

DeltaC_0 = 1.083e-03;

[dpb]

So have you tried rerunning nlinfit with those as initial guesses?

[Ronny Rives]

Yes, and I got the previous warning message. What means: Rank deficient, rank = 2, tol = 4.731182e-03? Where I can change rank or tolerance of nlinfit?

May be it is crucial to use some Robust option in the nlinfit? I tried before with robust option and two weight function ('andrews' and 'bisquare') but I got the same warning.

[dpb]

MOVED Ronny Answer to Comment...dpb

Finally I could get the estimated parameters without any warning. It was a scaling problem, one of the parameters were in the order of D = 1e-09 ~ 1e-10, while t_0 was approximately ~50. In addition, the two independent variables were in very different orders, time t = 5e + 04 and position z = 1e-03 ~ 1e-02. A simple change in the units of these variables avoided the previous numerical problems.

Thank you very much dpb, for all the help.

[dpb]

That scale difference was first thing that caught my eye when did get the data to look at, too. :)

That was the next thing I was going to suggest if simply a better starting point wasn't enough (and once I saw the data I didn't think likely). The scale difference shows up in the Jacobean=0 message and that wasn't dependent upon t almost certainly was the case of underflow for the exp(t) term--when t is 5E4 it takes a miniscule other prodect to keep the term from being identically zero with floating point precision.

Scaling of independent variables also can lead to the appearance of rank deficiency when comparing across columns -- the differences can be such that there is essentially no difference in numerical value because there isn't enough precision available to be able to store the difference--the small one just falls past the number of bits of precision in the larger so both columns appear to be almost identically the larger. There may be a small fraction of the elements towards the beginning of the time series that aren't but when get very far in the issue arises and overall the two look almost the same. The message doesn't say they are absolutely identical, only that they have a very high correlation.

That's why standardization is so often successful -- it takes the absolute value out and scales to something around +/-mean; the range being dependent upon the extremes from that mean, of course.

Or, as you did, simply re-scaling can work if it is just a magnitude issue like with the exponential term here...

Answer 2

[Ronny Rives]

Now I have another question. After applying the nlinfit algorithm to my experimental data, I obtained the estimated parameters of the model. Also, from CovB matrix I get the standard deviation (relative uncertainty) of these parameters. I would like to consider the uncertainty of my experimental data (which is known) in the uncertainty of my parameters. How could I do this?

Thanks in advance.