How do I fit a regression equation to find coefficients and exponents?

I'm currently working with a wide dataset, where I'm using a machine learning technique to select key identifiers (SIMPLS partial least squares). Then I want to use the identifiers and my outcome to create a predictive equation. I've tried a bunch of linear regression tools but they only find the predictor's coefficents, where I am trying to find the coefficients and the exponents. To get around this I'm trying to use 'nlinfit' to force the final equation into the desired form. This is where I'm having an issue, when I run the code I get the following error:

"The function you provided as the MODELFUN input has returned Inf or NaN values."

I've also tried inputting the model in the following form:

modelfun = 'y~(b1 + b2*x1.^b3 + b4*x2.^b5 + b6*x3.^b7)';

For reference my current data set is a 13x16 matrix, at it's largest it will be 24x35 matrix where the last column represents the outcome. Once the variable selection is complete (works without issue) the matrix is reduced to an nx4 matrix

Here is my code:

clear

clc

% Imports data and removed first text column
data = readtable('PMHS PLS Practice.xlsx',"textType","string");
data.SpecimenID = [];
% Splits data into independant and dependant 
% variables and normalizes values
X = data(:,1:15);
X = X{:,:};
Xnorm = normalize(X);
Y = data(:,16);
Y = table2array(Y);
Ynorm = normalize(Y);
% Performs Simpls PLS on normalized data returning
% the X scores and percent varience per variable
% for the first 5 latent variables
Cpart = cvpartition(13,"LeaveOut");
[~,~,SCR,~,~,PCTVAR,~,~] = plsregress(Xnorm,Ynorm,5,'cv',Cpart);
X_VAR = PCTVAR(1,:);
Y_VAR = PCTVAR(2,:);
% pareto(X_VAR)
% finds the dependant variable with the largest 
% contribution to the first 3 latent variables
[Var1,id1] = max(abs(SCR(:,1)));
[Var2,id2] = max(abs(SCR(:,2)));
[Var3,id3] = max(abs(SCR(:,3)));
% creates a matrix containing the selected variables 
X_reg = [data{:,id1} data{:,id2} data{:,id3}];
% Fits the data with a non-linear model with 
% initial coefficient guesses of beta0
modelfun = @(b,x) (b(1)+b(2)*x(:,1).^b(3)+b(4)*x(:,2).^b(5)+b(6)*x(:,3).^b(7));
beta0 = ones(1,7);
[coeff] = nlinfit(X_reg,Y,modelfun,beta0);

[Torsten]

We don't know your data.

But it is mandatory that x(:,1), x(:,2) and x(:,3) are vectors with only positive elements and without NaN or +/- Inf values.

Is this the case in your application ?

[Liam]

Yes, all of the vectors contain finite positive integers.

[Torsten]

Before calling nlinfit, try what you get when you call

Y_fit0 = modelfun(beta0,X_reg)

This is what nlinfit will do first and which results in the error message you get.

Does Y_fit0 contain Inf or NaN values as claimed by the error message ?

[Liam]

The code runs exactly the same with or without the the Y_fit0 line; all of the values in the Y_fit0 vector are approximately 1.3 if that makes a difference. I have copied all of the error messages below, for reference line 40 refers to the last line in the origional code. I believe that the first 3 errors are associated with what is happening in the curve fitting toolbox.

Error using nlinfit>checkFunVals

The function you provided as the MODELFUN input has returned Inf or NaN values.

Error in nlinfit>LMfit (line 596)

if funValCheck && ~isfinite(sse), checkFunVals(r); end

Error in nlinfit (line 284)

[beta,J,~,cause,fullr] = LMfit(X,yw, modelw,beta,options,verbose,maxiter);

Error in non_linear_model (line 40)

[coeff] = nlinfit(X_reg,Y,modelfun,beta0);

[Matt J]

It would help if you attached a .mat file containing Y and X_reg.

Also, are the exponent parameters b([3,5,7]) supposed to be non-negative? If so, it might be better to use a different fitting tool, one which allows you to specify lower bounds (e.g., lsqcurvefit or fit).

[Torsten]

Are you sure that there are no zero values or values < 0 in the vectors X_reg ?

Are you sure Y contains no NaN or +/- Inf values ?

If you are sure for both questions, then make a function

function y = modelfun(b,x)
  y = b(1)+b(2)*x(:,1).^b(3)+b(4)*x(:,2).^b(5)+b(6)*x(:,3).^b(7);
  if any(isnan(y)) || any(isinf(y))
    disp('Error in modelfun')
  end
end

instead of your function handle and check when and why the output vector y has Inf or NaN values in it.

[Liam]

@Matt J Here are .mat files containing the data for X_reg and Y. All of the values are positive integers, some relatively small but positive.

As to the exponents, I would assume them to be positive but they don't nessesarily need to be as the idea is to fit the relationship as close as possible. I'm not overly concerned with overfitting the data because the PLS is used to reduce the dimensionality prior to fitting the equation.

[Liam]

@Torsten using the function that you provided will only return the Y_fit0 values because it's unsupervised, ie I'm passing the function the initial guesses for the beta values and the x predictors, thus it will only return the initial properties for Y. See column vector below:

ans =

1.3500

1.3731

1.3516

1.3606

1.3627

1.3422

1.3692

1.3812

1.3784

1.3908

1.3571

1.3830

1.3913

[David Goodmanson]

Hi Liam, so x(:,1) etc. contain no zeros, is that correct?

[Liam]

Hi David, yes, all of the x(:,n) vectors are positive non-zero integers.

[Torsten]

using the function that you provided will only return the Y_fit0 values because it's unsupervised, ie I'm passing the function the initial guesses for the beta values and the x predictors, thus it will only return the initial properties for Y.

No. If you remove the line

modelfun = @(b,x) (b(1)+b(2)*x(:,1).^b(3)+b(4)*x(:,2).^b(5)+b(6)*x(:,3).^b(7));

in your code and add the function I wrote instead to your code, this function will be called every time by nlinfit if the solver wants to evaluate Yfit. And the error message will be displayed if any of the Yfit values is NaN or +/-Inf.

Maybe it helps further to print the b-vector in this case:

function y = modelfun(b,x)
y = b(1)+b(2)*x(:,1).^b(3)+b(4)*x(:,2).^b(5)+b(6)*x(:,3).^b(7);
if any(isnan(y)) || any(isinf(y))
    disp('Error in modelfun')
    b
end
end

[Liam]

I'm not sure that I exactly understand what you're saying? would I need to pass beta0 and X_reg to the function when I call nlinfit? like this?

[coeff] = nlinfit(X_reg,Y,modelfun(beta0,X_reg),beta0);

% Fits the data with a non-linear model with 
% initial coefficient guesses of beta0
beta0 = ones(1,7);
[coeff] = nlinfit(X_reg,Y,modelfun,beta0);

function y = modelfun(b,x)
y = b(1)+b(2)*x(:,1).^b(3)+b(4)*x(:,2).^b(5)+b(6)*x(:,3).^b(7);
if any(isnan(y)) || any(isinf(y))
    disp('Error in modelfun')
end
end

[Torsten]

[coeff] = nlinfit(X_reg,Y,@(b,x)modelfun(b,x),beta0);

instead of

[coeff] = nlinfit(X_reg,Y,modelfun,beta0);

[Matt J]

Here are .mat files containing the data for X_reg and Y.

Why two .mat files?

All of the values are positive integers, some relatively small but positive.

I see non-integer values in your attached files,

>> format long; Y(1), X_reg(1)

ans =

4.562586000000000e+03

ans =

0.062228000000000

In fact all X_reg are non-integers less than 1.

[Liam]

Okay so I tried that, and I recieved the same error as before, for some reason it just doesn't like the model function.

Error using nlinfit>checkFunVals

The function you provided as the MODELFUN input has returned Inf or NaN values.

Error in nlinfit>LMfit (line 596)

if funValCheck && ~isfinite(sse), checkFunVals(r); end

Error in nlinfit (line 284)

[beta,J,~,cause,fullr] = LMfit(X,yw, modelw,beta,options,verbose,maxiter);

Error in non_linear_model (line 40)

[coeff] = nlinfit(X_reg,Y,@modelfun,beta0);

[Torsten]

And what coefficients and what Yfit values did you get from "modelfun" at this stage ?

[Liam]

@Matt J

Why two .mat files?

I can't provide my raw data because it has identifiable donor health information. Since Y and X_reg have been preprocessed and filtered there isn't anything identifyible left. Also it's easier to provide the required info and eliminate the clutter.

I see non-integer values in your attached files'

You are correct, I was looking at the preprocessed data set but when I normalized it they became floating values.

Thanks for fminpleas I'll take a look at it!

[Liam]

@Torsten When applying the model as a function, the code doesn't execute the nlinfit line so the Yfit = Yfit0 and the coefficients remain zeros.

Yfit0 =

1.3500

1.3731

1.3516

1.3606

1.3627

1.3422

1.3692

1.3812

1.3784

1.3908

1.3571

1.3830

1.3913

[Matt J]

I can't provide my raw data because it has identifiable donor health information.

But why 2 .mat files instead of 1.

Answer 1

[Matt J]

fminspleas from the File Exchange (Download) fared better, but it looks like a highly ill-conditioned data set:

load('X_reg.mat')
load('Y.mat')

flist={1, @(b,x)x(:,1).^b(1), @(b,x)x(:,2).^b(2) , @(b,x)x(:,3).^b(3)};

[exps,coef]=fminspleas(flist,ones(1,3),X_reg,Y);

Warning: Rank deficient, rank = 3, tol = 1.539082e-12.

Warning: Rank deficient, rank = 3, tol = 9.749423e-13.

Warning: Rank deficient, rank = 3, tol = 2.738125e-13.

Warning: Rank deficient, rank = 3, tol = 4.562448e-13.

Warning: Rank deficient, rank = 3, tol = 2.633697e-13.

Warning: Rank deficient, rank = 3, tol = 3.101994e-13.

Warning: Rank deficient, rank = 3, tol = 3.745161e-13.

Warning: Rank deficient, rank = 3, tol = 3.900413e-13.

Warning: Rank deficient, rank = 3, tol = 2.828501e-13.

Warning: Rank deficient, rank = 3, tol = 3.704969e-13.

Warning: Rank deficient, rank = 3, tol = 2.850406e-13.

Warning: Rank deficient, rank = 3, tol = 3.122743e-13.

Warning: Rank deficient, rank = 3, tol = 2.781064e-13.

Warning: Rank deficient, rank = 3, tol = 2.855203e-13.

Warning: Rank deficient, rank = 3, tol = 2.706331e-13.

Warning: Rank deficient, rank = 3, tol = 2.717801e-13.

Warning: Rank deficient, rank = 3, tol = 2.854690e-13.

Warning: Rank deficient, rank = 3, tol = 2.848331e-13.

Warning: Rank deficient, rank = 3, tol = 2.704667e-13.

Warning: Rank deficient, rank = 3, tol = 2.409175e-13.

Warning: Rank deficient, rank = 3, tol = 2.643188e-13.

Warning: Rank deficient, rank = 3, tol = 2.525788e-13.

Warning: Rank deficient, rank = 3, tol = 2.539529e-13.

Warning: Rank deficient, rank = 3, tol = 2.792916e-13.

Warning: Rank deficient, rank = 3, tol = 2.631564e-13.

Warning: Rank deficient, rank = 3, tol = 2.728328e-13.

Warning: Rank deficient, rank = 3, tol = 2.672025e-13.

Warning: Rank deficient, rank = 3, tol = 2.548278e-13.

Warning: Rank deficient, rank = 3, tol = 2.646537e-13.

Warning: Rank deficient, rank = 3, tol = 2.597743e-13.

Warning: Rank deficient, rank = 3, tol = 2.603938e-13.

Warning: Rank deficient, rank = 3, tol = 2.636856e-13.

Warning: Rank deficient, rank = 3, tol = 2.613920e-13.

Warning: Rank deficient, rank = 3, tol = 2.626564e-13.

Warning: Rank deficient, rank = 3, tol = 2.617220e-13.

Warning: Rank deficient, rank = 3, tol = 2.627517e-13.

Warning: Rank deficient, rank = 3, tol = 2.623634e-13.

Warning: Rank deficient, rank = 3, tol = 2.636579e-13.

Warning: Rank deficient, rank = 3, tol = 2.626151e-13.

Warning: Rank deficient, rank = 3, tol = 2.630883e-13.

Warning: Rank deficient, rank = 3, tol = 2.626463e-13.

Warning: Rank deficient, rank = 3, tol = 2.628092e-13.

Warning: Rank deficient, rank = 3, tol = 2.626155e-13.

Warning: Rank deficient, rank = 3, tol = 2.629409e-13.

Warning: Rank deficient, rank = 3, tol = 2.628521e-13.

Warning: Rank deficient, rank = 3, tol = 2.626943e-13.

Warning: Rank deficient, rank = 3, tol = 2.622067e-13.

Warning: Rank deficient, rank = 3, tol = 2.627968e-13.

Warning: Rank deficient, rank = 3, tol = 2.625353e-13.

Warning: Rank deficient, rank = 3, tol = 2.625002e-13.

Warning: Rank deficient, rank = 3, tol = 2.624827e-13.

Warning: Rank deficient, rank = 3, tol = 2.636432e-13.

Warning: Rank deficient, rank = 3, tol = 2.631713e-13.

Warning: Rank deficient, rank = 3, tol = 4.432535e-13.

Warning: Rank deficient, rank = 3, tol = 4.432535e-13.

Warning: Rank deficient, rank = 3, tol = 4.164941e-13.

Warning: Rank deficient, rank = 3, tol = 5.000336e-13.

Warning: Rank deficient, rank = 3, tol = 4.987812e-13.

Warning: Rank deficient, rank = 3, tol = 4.468363e-13.

Warning: Rank deficient, rank = 3, tol = 5.283489e-13.

Warning: Rank deficient, rank = 3, tol = 4.646775e-13.

Warning: Rank deficient, rank = 3, tol = 5.457478e-13.

Warning: Rank deficient, rank = 3, tol = 5.054967e-13.

Warning: Rank deficient, rank = 3, tol = 5.469517e-13.

Warning: Rank deficient, rank = 3, tol = 5.469517e-13.

Warning: Rank deficient, rank = 3, tol = 5.653845e-13.

Warning: Rank deficient, rank = 3, tol = 5.420883e-13.

Warning: Rank deficient, rank = 3, tol = 5.196659e-13.

Warning: Rank deficient, rank = 3, tol = 5.498905e-13.

Warning: Rank deficient, rank = 3, tol = 5.498905e-13.

Warning: Rank deficient, rank = 3, tol = 5.255071e-13.

Warning: Rank deficient, rank = 3, tol = 5.494949e-13.

Warning: Rank deficient, rank = 3, tol = 5.379439e-13.

Warning: Rank deficient, rank = 3, tol = 5.386419e-13.

Warning: Rank deficient, rank = 3, tol = 5.328367e-13.

Warning: Rank deficient, rank = 3, tol = 5.427387e-13.

Warning: Rank deficient, rank = 3, tol = 5.372832e-13.

Warning: Rank deficient, rank = 3, tol = 5.430067e-13.

Warning: Rank deficient, rank = 3, tol = 5.389417e-13.

Warning: Rank deficient, rank = 3, tol = 5.387856e-13.

Warning: Rank deficient, rank = 3, tol = 5.411788e-13.

Warning: Rank deficient, rank = 3, tol = 5.393915e-13.

Warning: Rank deficient, rank = 3, tol = 5.402901e-13.

Warning: Rank deficient, rank = 3, tol = 5.395469e-13.

Warning: Rank deficient, rank = 3, tol = 5.398723e-13.

Warning: Rank deficient, rank = 3, tol = 5.398540e-13.

Warning: Rank deficient, rank = 3, tol = 5.399638e-13.

Warning: Rank deficient, rank = 3, tol = 5.397106e-13.

Warning: Rank deficient, rank = 3, tol = 5.406220e-13.

Warning: Rank deficient, rank = 3, tol = 5.405250e-13.

Warning: Rank deficient, rank = 3, tol = 5.392107e-13.

Warning: Rank deficient, rank = 3, tol = 5.400403e-13.

Warning: Rank deficient, rank = 3, tol = 5.395216e-13.

Warning: Rank deficient, rank = 3, tol = 5.398500e-13.

Warning: Rank deficient, rank = 3, tol = 5.395714e-13.

Warning: Rank deficient, rank = 3, tol = 5.397015e-13.

Warning: Rank deficient, rank = 3, tol = 5.397013e-13.

Warning: Rank deficient, rank = 3, tol = 5.390639e-13.

Warning: Rank deficient, rank = 3, tol = 5.393260e-13.

Warning: Rank deficient, rank = 3, tol = 5.393357e-13.

Warning: Rank deficient, rank = 3, tol = 5.394242e-13.

Warning: Rank deficient, rank = 3, tol = 5.389746e-13.

Warning: Rank deficient, rank = 3, tol = 5.392055e-13.

Warning: Rank deficient, rank = 3, tol = 5.391447e-13.

Warning: Rank deficient, rank = 3, tol = 5.391693e-13.

Warning: Rank deficient, rank = 3, tol = 5.391516e-13.

Warning: Rank deficient, rank = 3, tol = 5.390076e-13.

Warning: Rank deficient, rank = 3, tol = 5.390934e-13.

Warning: Rank deficient, rank = 3, tol = 5.390303e-13.

Warning: Rank deficient, rank = 3, tol = 5.390635e-13.

Warning: Rank deficient, rank = 3, tol = 5.392366e-13.

Warning: Rank deficient, rank = 3, tol = 5.390566e-13.

Warning: Rank deficient, rank = 3, tol = 5.391245e-13.

Warning: Rank deficient, rank = 3, tol = 5.390133e-13.

Warning: Rank deficient, rank = 3, tol = 5.390656e-13.

Warning: Rank deficient, rank = 3, tol = 5.390554e-13.

Warning: Rank deficient, rank = 3, tol = 5.389963e-13.

[b1,b2,b4,b6]=dealThem(coef)

b1 = 1.0624e+04

b2 = 1.7367e+06

b4 = -116.6191

b6 = -3.0221e+15

[b3,b5,b7]=dealThem(exps)

b3 = 2.6098

b5 = -2.8463

b7 = 9.4736

function varargout=dealThem(z)
   varargout=num2cell(z);
end

[Liam]

So I've had a little time to play around with fminspleas and it definitely works for what I'm trying to do, that being said I'm going to look around and see if theres anything that handles highly correlated data a little better.

Weirdly enough, I tried a few different things on fminspleas (weighted and unweighted) with the same data that I provided and I always got a different answer than you. Though yours seems to achieve a better tolerance.

Warning: Rank deficient, rank = 3, tol = 2.020925e-07.

Warning: Rank deficient, rank = 3, tol = 6.902514e-08.

Warning: Rank deficient, rank = 3, tol = 4.270030e-07.

Warning: Rank deficient, rank = 2, tol = 2.980757e+00.

Warning: Rank deficient, rank = 3, tol = 4.465878e-03.

Warning: Rank deficient, rank = 2, tol = 2.938562e+00.

Warning: Rank deficient, rank = 3, tol = 3.779822e-02.

Warning: Rank deficient, rank = 3, tol = 3.529404e-02.

Warning: Rank deficient, rank = 3, tol = 1.624543e-02.

Warning: Rank deficient, rank = 3, tol = 7.019358e-02.

Warning: Rank deficient, rank = 3, tol = 5.739297e-04.

Warning: Rank deficient, rank = 3, tol = 3.844059e-04.

Warning: Rank deficient, rank = 3, tol = 8.231624e-05.

Warning: Rank deficient, rank = 3, tol = 2.292359e-04.

Warning: Rank deficient, rank = 3, tol = 2.292359e-04.

Warning: Rank deficient, rank = 3, tol = 7.355443e-04.

Warning: Rank deficient, rank = 3, tol = 5.833903e-05.

Warning: Rank deficient, rank = 3, tol = 4.428108e-05.

Warning: Rank deficient, rank = 3, tol = 8.290589e-05.

Warning: Rank deficient, rank = 3, tol = 1.123903e-04.

Warning: Rank deficient, rank = 3, tol = 6.118206e-05.

Warning: Rank deficient, rank = 3, tol = 5.337120e-05.

Warning: Rank deficient, rank = 3, tol = 5.616255e-05.

Warning: Rank deficient, rank = 3, tol = 5.616255e-05.

Warning: Rank deficient, rank = 3, tol = 5.814098e-05.

Warning: Rank deficient, rank = 3, tol = 3.288430e-05.

Warning: Rank deficient, rank = 3, tol = 2.534922e-05.

Warning: Rank deficient, rank = 3, tol = 2.868990e-05.

Warning: Rank deficient, rank = 3, tol = 3.289206e-05.

Warning: Rank deficient, rank = 3, tol = 3.410611e-05.

Warning: Rank deficient, rank = 3, tol = 3.212528e-05.

Warning: Rank deficient, rank = 3, tol = 3.507096e-05.

Warning: Rank deficient, rank = 3, tol = 2.276614e-05.

Warning: Rank deficient, rank = 3, tol = 2.989308e-05.

Warning: Rank deficient, rank = 3, tol = 2.207071e-05.

Warning: Rank deficient, rank = 3, tol = 2.651110e-05.

Warning: Rank deficient, rank = 3, tol = 2.348670e-05.

Warning: Rank deficient, rank = 3, tol = 2.522671e-05.

Warning: Rank deficient, rank = 3, tol = 2.048854e-05.

Warning: Rank deficient, rank = 3, tol = 2.009223e-05.

Warning: Rank deficient, rank = 3, tol = 2.001502e-05.

Warning: Rank deficient, rank = 3, tol = 1.720824e-05.

Warning: Rank deficient, rank = 3, tol = 2.155504e-05.

Warning: Rank deficient, rank = 3, tol = 1.851085e-05.

Warning: Rank deficient, rank = 3, tol = 1.913940e-05.

Warning: Rank deficient, rank = 3, tol = 1.791282e-05.

Warning: Rank deficient, rank = 3, tol = 1.922028e-05.

Warning: Rank deficient, rank = 3, tol = 1.912505e-05.

Warning: Rank deficient, rank = 3, tol = 1.807495e-05.

Warning: Rank deficient, rank = 3, tol = 1.843179e-05.

Warning: Rank deficient, rank = 3, tol = 1.973779e-05.

Warning: Rank deficient, rank = 3, tol = 1.913021e-05.

Warning: Rank deficient, rank = 3, tol = 1.875068e-05.

Warning: Rank deficient, rank = 3, tol = 1.875068e-05.

Warning: Rank deficient, rank = 3, tol = 1.960611e-05.

Warning: Rank deficient, rank = 3, tol = 1.940055e-05.

Warning: Rank deficient, rank = 3, tol = 1.940055e-05.

Warning: Rank deficient, rank = 3, tol = 1.917227e-05.

Warning: Rank deficient, rank = 3, tol = 1.934204e-05.

Warning: Rank deficient, rank = 3, tol = 1.892093e-05.

Warning: Rank deficient, rank = 3, tol = 1.923863e-05.

Warning: Rank deficient, rank = 3, tol = 1.909304e-05.

Warning: Rank deficient, rank = 3, tol = 1.904800e-05.

Warning: Rank deficient, rank = 3, tol = 1.920230e-05.

Warning: Rank deficient, rank = 3, tol = 1.909907e-05.

Warning: Rank deficient, rank = 3, tol = 1.915905e-05.

Warning: Rank deficient, rank = 3, tol = 1.911747e-05.

Warning: Rank deficient, rank = 3, tol = 1.915911e-05.

Warning: Rank deficient, rank = 3, tol = 1.914215e-05.

Warning: Rank deficient, rank = 3, tol = 1.921714e-05.

Warning: Rank deficient, rank = 3, tol = 1.915963e-05.

Warning: Rank deficient, rank = 3, tol = 1.918753e-05.

Warning: Rank deficient, rank = 3, tol = 1.916345e-05.

Warning: Rank deficient, rank = 3, tol = 1.917339e-05.

Warning: Rank deficient, rank = 3, tol = 1.916301e-05.

Warning: Rank deficient, rank = 3, tol = 1.916336e-05.

Warning: Rank deficient, rank = 3, tol = 1.915697e-05.

Warning: Rank deficient, rank = 3, tol = 1.914124e-05.

Warning: Rank deficient, rank = 3, tol = 1.916045e-05.

Warning: Rank deficient, rank = 3, tol = 1.915290e-05.

Warning: Rank deficient, rank = 3, tol = 1.916424e-05.

Warning: Rank deficient, rank = 3, tol = 1.916424e-05.

Warning: Rank deficient, rank = 3, tol = 1.915764e-05.

Warning: Rank deficient, rank = 3, tol = 1.916416e-05.

coef

b1 = 6.6081e+03

b2 = 1.7901e+07

b4 = -0.0929

b6 = -8.1415e-08

exps

b3 = 3.8416

b5 = -6.0307

b7 = -6.5186

[Matt J]

Weirdly enough, I tried a few different things on fminspleas (weighted and unweighted) with the same data that I provided and I always got a different answer than you.

Not weird at all. It's to be expected, based on my comment to Alex below.

[Matt J]

If you pre-normalize the columns of X_reg, fminspleas gives the same results as Alex and my conditioning test shows a much better condition number on the solution for the coefficients:

load('X_reg.mat')
load('Y.mat')

flist={1, @(b,x)x(:,1).^b(1), @(b,x)x(:,2).^b(2) , @(b,x)x(:,3).^b(3)};

s=max(X_reg);
n=size(Y,1);

[exponents,coefficients]=fminspleas(flist,ones(1,3),X_reg./s,Y);

coefficients(2:end)=coefficients(2:end)./(s.^exponents)';

format longG
exponents,coefficients

exponents = 1×3

109.691602052618 -3.07105187581003 12.1425625465001

coefficients = 4×1

1.0e+00 * 12988.1091748898 1.53571603411212e+126 -96.9999246561493 -9.29878054932722e+18

A=(X_reg./s).^exponents; A=[ones(n,1),A];

cond(A)

ans = 
          24.9901708703067

Answer 2

[Alex Sha]

Although the results may seem strange, mathematically speaking, the result below is the best one:

Sum Squared Error (SSE): 875229.002284955
Root of Mean Square Error (RMSE): 259.47120816783
Correlation Coef. (R): 0.945217846153423
R-Square: 0.893436776686916

Parameter	Best Estimate      	Std. Deviation       	Confidence Bounds[95%] 
---------	-------------      	--------------       	--------------------------------
b1       	12988.1117515585   	6.29998367140824     	[12972.6962468509, 13003.5272562661]
b2       	1.533510681096E126 	0.070612379218076    	[1.533510681096E126, 1.533510681096E126]
b3       	109.691046465045   	0.167406667907439    	[109.281417105381, 110.100675824708]
b4       	-97.0000950427331  	1048.74886499901     	[-2663.19612168365, 2469.19593159819]
b5       	-3.07105087006786  	14192.7795542942     	[-34731.5515429606, 34725.4094412205]
b6       	-9.2988265393269E18	1.40462270054153E-15 	[-9.2988265393269E18, -9.2988265393269E18]
b7       	12.1425640359073   	3.30305643552503E-124	[12.1425640359073, 12.1425640359073]

[Matt J]

the result below is the best one

One of the best ones perhaps. I don't think there is a unique best fit.

[Alex Sha]

Should be unique, if ignoring the error caused by the accuracy of numerical computation

Sum Squared Error (SSE): 875229.002279992
Root of Mean Square Error (RMSE): 259.471208167095
Correlation Coef. (R): 0.945217846153723
R-Square: 0.893436776687483

Parameter	Best Estimate       
---------	-------------       
b1       	12988.1075276193    
b2       	1.53583923988692E126
b3       	109.691633072097    
b4       	-96.9998016037511   
b5       	-3.07105262671452   
b6       	-9.2988221706396E18 
b7       	12.142564075687 

[Matt J]

Should be unique, if ignoring the error caused by the accuracy of numerical computation

I don't know why you think so. Even if we assume that your exponents b3,b5,b7 are correct and true, the condition number on the remaining unknowns is terrible:

load X_reg;

b3=109.691633072097;
b5=-3.07105262671452;
b7=12.142564075687;

A=X_reg.^[b3,b5,b7]; A(:,end+1)=1;

cond(A)

ans = 4.2252e+125

[Alex Sha]

Hi, just discussion!

If there are various solutions as you said, what are the objective function values (SSE) Corresponding to those solutions? equal to or even less than the result I shown above, that is: SSE=875229.002279992, if your results are bad than this, it menas the results you get are all local solutions, not global one.

There is only one global optimization solution, like I provied previously.

Checking one result you given:

b1 = 1.0624e+04
b2 = 1.7367e+06
b4 = -116.6191
b6 = -3.0221e+15
b3 = 2.6098
b5 = -2.8463
b7 = 9.4736

The SSE value is: 1432079.63835075, it is much higher (poor) than 875229.002279992.

The reason why you get different solution each times computation is maybe the optimization algorithms you uesd are local optimization algorithms (LOA), not global optimization algorithms (GOA), LOA has no ability, and does not guarantee a global optimal solution.

Anyway, if you can get SSE value near to 875229.002279992, you will find that all parameter values will be fixed.

[Matt J]

If there are various solutions as you said, what are the objective function values (SSE)

Probably very similar to what you got. My solution may be local and your solution may be global, but that does not mean the global solution is unique.

[Alex Sha]

If possible, would you please be kind to show me one more global solution other than mine? so I can make some comparisons and find out why.