How to tweak an equation to properly fit a dataset using lsqcurvefit?

I have an equation that I'm working on that has 3 independent variables and may need to include additional fractional terms to properly fit a dataset. Here is the original equation:

Units:

y in m/s

and in m

in Pa

in

is unitless (Poisson's ratio)

I've tried setting the parameters to be used in fitting as such:

where and are the parameters.

However, this setup was not sufficient for the equation to be fitted to my dataset, as there were huge differences.

My question is, how can I know which fractional terms to add to this equation to be used in the fitting process, while also conserving the units of the output? For instance, I tried modifying the equation like this:

but the fit did not change, as far as I can tell visually.

I am using lsqcurvefit for fitting, here is my code:

independentVariables = [x1 x2 x3];

y = @(p, x) p(1) .* sqrt((2 .* x(:,1)) + p(3) .* (x(:,2) ./ x(:,3)) .^ p(4) ./ ...
                           (c1 .* (3 + c2) .* (1 - ((x(:,2) - x(:,3)) ./ x(:,2)) .^ p(2))));

% Initial guess for parameters, lower and upper bounds
initialGuess = [1 1 1];
lb = [];
ub = [];

% Call lsqcurvefit
parameters_y = lsqcurvefit(y, initialGuess, independentVariables, ydata(:,1),lb,ub);

Here is my dataset, I also added it as an attachment to this post:

ydata =

36165.2399652261
36832.2283727431
38244.8839743203
39532.2710645211
38985.1533416210
34225.9039660907
30164.7430624110
26359.0764634342
22772.0684138657
24623.4765316140
8877.15808522203
8717.70659273395
8558.63955174474
8603.93260894558
8659.68585200240
6188.18034995559
6123.01438079126
6057.56501843796
5991.56509717007
5924.71769692532
4872.24385495395
4852.85576311786
4832.76067401530
4811.90659475242
4790.23277171101

ydata = reshape(ydata,5,[]) % ydata is actually a 2D matrix where rows vary with x2 and columns vary with x3.

And finally here is how my fit looks right now:

Different angle:

Any help would be appreciated.

Edit: Added the independent variables, where columns in independentVariables are .

Edit 2: Added units of c1, c2, x1, x2, x3 and y.

Edit 3: Corrected the variables in the post to be in line with the .mat files.

[dpb]

load ydata
load independentVariables
whos

  Name                       Size            Bytes  Class     Attributes

  ans                        1x49               98  char                
  cmdout                     1x33               66  char                
  independentVariables      25x3               600  double              
  ydata                     25x1               200  double              

Where are x1, x2, x3? Nothing anybody can do without them.

plot(ydata)

Just an observation, the plateau with a sudden break such as shown is not going to be fit well with a continuous function; that looks like a discrete parameter model may be in order at whatever it is that differentiates the first and second halves of the observations.

Also, with three independent variables, having only 25 observations isn't going to be much data in each.

[Juneight]

Apologies, attached a new variable called independentVariables to the post.

I don't know if it's a good practice or not, but I stored the results as a vector rather than a 2D matrix in ydata. That is what accounts for those sudden changes from 10th to 11th point, for example.

I can't add anymore genuine observations, but would you say interpolating to have more datapoints in ydata be helpful during fitting?

[William Rose]

@Juneight,

Good comments from @dpb.

You asked "how can I know which fractional terms to add to this equation to be used in the fitting process, while also conserving the units of the output?" Please specify the units for y, x1, x2, x3 if you want advice on consistency of units. Your equation includes the constannts 2, 3, c1, c2. Do those constants have units attached?

What are the values of c1, c2? If c1, c2 were obtained by a previous attempt to fit the data, then you'd want to fit them at the same time as p(). As @dpb noted fitting four parameters to 25 points is a stretch, and six is definitely too many, I think.

When you try different versions of the equation, to see what works, you are, in a way, fitting even more parameters, in the many-dimensional space of "equations that could work". The risk of overfitting increases when you allow youself the freedom to try various equations to fit your data.

[Juneight]

Hello @William Rose,

I added the units for the variables to the post.

c1 and c2's values were obtained previously without fitting. They are the same for all 25 observations. I thought their exact values would be irrelevant for this question but I can add them tomorrow if it's necessary.

What would you recommend for the number of parameters to be used in such an equation? I have limited experience in fitting, previously in this project, I fitted a different equation with 3 independent variables and 3 parameters with about 5% error using lsqcurvefit.

[dpb]

load ydata
load independentVariables
whos

  Name                       Size            Bytes  Class     Attributes

  ans                        1x49               98  char                
  cmdout                     1x33               66  char                
  independentVariables      25x3               600  double              
  ydata                     25x1               200  double              

format shortg
v=independentVariables;
[v (v(:,2)-v(:,3))./v(:,2) ydata]

ans = 25×5

1.0e+00 * 0.0125 0.001 5.4149e+07 -5.4149e+10 36165 0.0125 0.00325 1.7123e+08 -5.2686e+10 36832 0.0125 0.0055 2.892e+08 -5.2582e+10 38245 0.0125 0.00775 3.9467e+08 -5.0925e+10 39532 0.0125 0.01 4.3402e+08 -4.3402e+10 38985 0.0325 0.001 2.4845e+07 -2.4845e+10 34226 0.0325 0.00325 7.6139e+07 -2.3428e+10 30165 0.0325 0.0055 1.2079e+08 -2.1962e+10 26359 0.0325 0.00775 1.5838e+08 -2.0436e+10 22772 0.0325 0.01 2.0318e+08 -2.0318e+10 24623

Wow!!! With over 10 orders of magnitude between x1,x2 and x3, the ratio term is going to saturate the model response; x1 in the numerator will essentialy turn the very large positive exponent into a very large negative one, but the relative size remains and 10^10 squared is 10^20 and everything else pales in comparison.

For such a model, one would presume there must be some physical basis behind choosing such, where does this form derive? 1E8 m --> 62,000 miles as compared to mils for x1, x2. What in the world are we measuring here?

[William Rose]

@Juneight,

[edit: fix typo]

@dpb is very smart and observant. He looked at the independent variables and immediately noticed that and , by a factor of 10^10 or so. Your models are

(model 1)

or

(model 2)

Both models include the term , and @dbp recognizes that this is not significantly different from . Do you have a physical reason to believe that (x2-x3) is a factor in your model? It seems implausible.

Since is of order 10^10, is very sensitive to p(2), in a non-linear way. If p(2)>0, then .

If, in addition, p(4)<0 and abs(p(3)) is not tiny, then

,

and now model 2 is, to a very good approximation,

The model above is unfittable (ill-conditioned), because p(1) and p(3) cannot be independently determined (a change in one can be perfectly offset by a change in the other), and likewise for p(2) and p(4).

In your code, you did not specify lower and upper bounds for the parameters. It is a good idea to include bounds for all parameters being fitted: non-negative, less than a million, or something. This is not hard, if the parameters have some physical meaning.

[William Rose]

@Juneight,

You asked, "What would you recommend for the number of parameters to be used in such an equation?".

That is a very good question with no simple answer. There is no specific number of parameters that are recommended. I have read (I forget where) that one should have at least 10 data points per fitted parameter. One could question the justification for this recommendation. I think it is not a bad place to start. Of course you want the data points to span ranges that are significantly larger than the measurement noise in their respective dimensions. A good approach, which takes some effort, is to compute error bars for the parameter estimates. If the 95% error bar for parameter P spans zero, then you probably should not keep parameter P. There are various approaches to estimating whether added parameters are justified in a model: Akaike information criterion, Bayesian information criterion, Rissanen's minimum description length are well known.

[Juneight]

Sorry @dpb and @William Rose, I made a mistake while preparing the question post and changed the naming of the independent variables, they should make more sense now, I also updated the post but the independent variable that was previously named should have been , as its in Pascals and it is the third column in the independentVariables.

"In your code, you did not specify lower and upper bounds for the parameters. It is a good idea to include bounds for all parameters being fitted: non-negative, less than a million, or something. This is not hard, if the parameters have some physical meaning."

Noted, I will add them.

Answer 1

[William Rose]

@Juneight,

{Edit:

Dimensional analysis needs changing because now you say x3 is [Pa] and before you said x3 is [m].

You said "the independent variable that was previously named x1 should have been x3, as its in Pascals and it is the third column in the independentVariables". Are you proposing that we exchange x1 and x3 in the model equation? Then my previous model analysis is now incorrect (and was a waste of time), and my model equations are incorrect. Arrgh!

When I compare youre new model 1 to the old model 1, I see that you did not just exchage x1 and x3 in the equation. What was x1 is now x3. What was x2 is now x1. What was x3 is now x2. @Juneight, am I going crazy? I'm definitely getting grumpy. Please clarify your comments regarding the exchanges of variables.}

I would not use interpolation to add points. This will not be a relaible way to improve the fit, since there seem to be non-linearities in the relationship, which we cannot account for in the interpolation. Also, this would be interpolation in 3 dimensions of x, which introduces more chance for error than interpolation in 1 or 2 dimensions.

The units for p(2) in model 1, and for p(2) and p(4) in model 2, are dimensionless. p(3) in model 2 has dimensions of meters. p(1) in models 1 and 2 is dimensionless. What is the physical interpretation of p(1)? It is nice when model parameters have physical meaning.

load('ydata'); 
load('independentVariables');
xdata=independentVariables;
% Make plots showing distribution of the independent variables
figure; subplot(131); plot(xdata(:,1),xdata(:,2),'r*'); xlabel('X1'); ylabel('X2'); grid on
subplot(132); plot(xdata(:,2),xdata(:,3),'g*'); xlabel('X2'); ylabel('X3'); grid on
subplot(133); plot(xdata(:,3),xdata(:,1),'b*'); xlabel('X3'); ylabel('X1'); grid on

Let's try fitting y=f(x1,x2,x3) with the model 1 equation, which is

[Insert the new equation which @Juneight substituted.]

In your first model (before model 1), you had 2 in place of p(2).

{The next several lines may not be corerct since @Juneight has redefined the variables x1, x2, x3.}

Therefore I assume that p(2) is probably a positive number, but maybe not 2. But since x3>>x2, and since x1 is positive, and c1 and c2 are positive, you will get a negative value under the square root, i.e. you'll get a complex number from this model. Therefore I suspect there is something wrong with your x3 or x2 data. Maybe x2 is in units of meters, and x3 is in nanometers. Or the model is wrong.

c1=1e3;                     % density (kg/m^3)
c2=0.33;                    % Poisson's ratio (non-dim)
%y = @(p, x) p(1)*sqrt(2*x(:,1)./(c1*(3+c2)*(1-((x(:,2)-x(:,3))./x(:,2)).^p(2))));
% What was x1 is now x3.  What was x2 is now x1.  What was x3 is now x2.
y = @(p, x) p(1)*sqrt(2*x(:,3)./(c1*(3+c2)*(1-((x(:,1)-x(:,2))./x(:,1)).^p(2))));
p0 = [1 2];                 % initial guess for p()
lb = [1e-2 0.3];            % lower bounds for p()
ub = [1e6 4];               % upper bounds for p()
% Estimate p and get the sum squared error (SSE)
[pEst,SSE] = lsqcurvefit(y, p0, xdata, ydata, lb, ub);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

fprintf('p1=%.3f, p2=%.3f, SSE=%.2e.\n',pEst(1),pEst(2),SSE)

p1=47.782, p2=4.000, SSE=3.74e+09.

[dpb]

Fixing the typo that removed the exponent term...

load('ydata'); 
load('independentVariables');
xdata=independentVariables;
c1=1e3;                     % density (kg/m^3)
c2=0.33;                    % Poisson's ratio (non-dim)
y = @(p, x) p(1)*sqrt(2*x(:,3)./(c1*(3+c2)*(1-((x(:,1)-x(:,2))./x(:,1)).^p(2))));
p0 = [1 2];                 % initial guess for p()
lb = [1e-2 0.3];            % lower bounds for p()
ub = [1e6 4];               % upper bounds for p()
% Estimate p and get the sum squared error (SSE)
[pEst,SSE] = lsqcurvefit(y, p0, xdata, ydata, lb, ub);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

fprintf('p1=%.3f, p2=%.3f, SSE=%.2e.\n',pEst(1),pEst(2),SSE)

p1=47.782, p2=4.000, SSE=3.74e+09.

Fixing up the model to match the revised definition of variables does then find a solution, but the SSE is a big number, so the fit probably isn't all that good. Am too lazy to try to figure out how to create the surface plot from the one-vector so will just plot response and error versus observation and see what that most rudimentary of diagnostics shows...

x=[1:numel(ydata)];
yfit=y(pEst,xdata);
res=ydata-yfit;
plot(x,[yfit(:) ydata(:) res(:)],'*')
legend('Fit','Data','Residual')

shows the fit magnitude is at least in the neighborhood. Now the model can at least begin to handle the disparity between the first/last half of the observations although there's still not nearly enough of a discontinuity introduced for those first points and the slope between the 5th to 10th points is reversed from the actual data trend so something is definitely not being accounted for there.

It is interesting it did come up with 2.000 for p(2); not sure what (if anything) that signifies. Well @William Rose had sharp-enough eyes to miss the fixup I missed in correcting to the revised location of the independent variables explains that anomaly. Now, the exponent is again bounded; we can try to expand its upper bound.

figure
plot(x,[yfit(:) ydata(:)],'*')
legend('Fit','Data')

Plotting without the residuals that throw up big negative numbers shows more clearly the model is a series of roughly parallel lines with a very slight slope as the observation number increases that mimics the last three sections moderately well, but again the first/second halves are totally disparate in behavior it appears and the slope of the fit even for the subsections is backwards from the observed data trend.

ub = [1e6 1E2];               % upper bounds for p()
% Estimate p and get the sum squared error (SSE)
[pEst,SSE] = lsqcurvefit(y, p0, xdata, ydata, lb, ub);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

fprintf('p1=%.3f, p2=%.3f, SSE=%.2e.\n',pEst(1),pEst(2),SSE)

p1=74.559, p2=49.592, SSE=2.45e+09.

yfit=y(pEst,xdata);
res=ydata-yfit;
figure
plot(x,[yfit(:) ydata(:)],'*')
legend('Fit','Data')

And now, we're back to an earlier result in which p(2) ~ 50 which gives the exponential term way too much weight and, with the constraint that C1, C2 >0 all those slopes are positive...now the very first set of points are, but those thereafter all all negative.

Would need to go to the effort to now work on the residuals plot against actual variables to see how those behave and, it got me too confused about how to rearrange the 25x3 array to get a gridded output that I didn't have time to wrap head around to get to those surface plots...

The net conclusion to me is the model is simply not representative of the data trends -- either the model doesn't describe the actual measurements from which it came or, perhaps there are one or more confounding variables not measured and hence not included in the model that cause the observed behavior.

OTOMH, I don't have a klew as to what one could introduce into anything at all resembling the posited model and having absolutely no context from which to even begin to derive some physical model leaves one grasping at air...

[William Rose]

@dpb, you are most patient.

@dbp, I think your reviised model equation has an typo:

y = @(p, x) p(1)*sqrt(2*x(:,3)./(c1*(3+c2)*(1-((x(:,1)-x(:,1))./x(:,1)).^p(2))));

should be

y = @(p, x) p(1)*sqrt(2*x(:,3)./(c1*(3+c2)*(1-((x(:,1)-x(:,2))./x(:,1)).^p(2))));

When I run the model with that correction, I get a different solution, not surprisingly. And the solution has p2 up against its upper limit of 4. Maybe the upper limit needs revising. I chose 4 as the upper limit for p2 because in "model 0", there were not fitted parameters, and p2 (which is a power) was just 2.

[William Rose]

@Juneight,

I just noticed that the model of @dpb, which has what I suspect is a typo, fits better (has a lower SSE) than my model, with the "typo" corrected. His model is

y = @(p, x) p(1)*sqrt(2*x(:,3)./(c1*(3+c2)*(1-((x(:,1)-x(:,1))./x(:,1)).^p(2))));

and mine is

y = @(p, x) p(1)*sqrt(2*x(:,3)./(c1*(3+c2)*(1-((x(:,1)-x(:,2))./x(:,1)).^p(2))));

If you look at his model, it simplifies to

y = @(p, x) p(1)*sqrt(2*x(:,3)./(c1*(3+c2)*(1-0^p(2))));

so p(2) does nothing, and x(1) and x(2) play no role, and the model simplifies further to

y = @(p, x) p(1)*sqrt(2*x(:,3)./(c1*(3+c2)));

Now we know why the fitted p(2) in @dpb's model was 2, which was the initial guess: because it plays no role, so it never changed. The simple model gave a lower SSE than the model involving p(2) and x(1) and x(2). Maybe we chose a bad intial guess for the intial model.

[dpb]

"Maybe we chose a bad intial guess for the intial model."

Think??? :)

Clearly the model doesn't represent these data -- maybe it fits some specific situation for which it was developed, but something is far different in the data set at hand and that model.

But, as noted before, we have insufficient knowledge of anything and everything to have any hope at all of understanding what it is, specifically are looking at and what any of it means (if anything) and how the data were collected, etc., etc., etc., ...

One thing obvious from your prior plots of the range of the indepenent variables is that the P term is quite constrained for most of the other variables; only one dataset traverses the whole range. That leaves a large area of that plane uncovered for estimation. Now, it may very well be that the other pressures can't be obtained for those other conditions for physical reasons...

[dpb]

surf

load ydata
load independentVariables
y=ydata; clear ydata
x=independentVariables; clear independentVariables
x1=x(1:5:end,1); x2=x(1:5,2);
subplot(1,3,1), plot(x(:,1),x(:,2),'*'), xlabel('x1'), ylabel('x2')
subplot(1,3,2), plot(x(:,1),x(:,3),'*'), xlabel('x1'), ylabel('x3')
subplot(1,3,3), plot(x(:,2),x(:,3),'*'), xlabel('x2'), ylabel('x3')

figure
z=reshape(y,5,5);
surf(x1,x2,z), xlabel('x1'), ylabel('x2')

Have to leave unfinished -- x1, x2 are uniform grids, but x3 is not, can't use surf() but will have to use scattered interpolant because there isn't a uniform mesh.

numel(unique(x(:,3)))

ans = 25

In fact, every pressure point is unique...

[William Rose]

“Think??” Makes me chuckle :) Interesting plot and interesting that all pressures are unique.