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

26 次查看(过去 30 天)
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);
  19 个评论
Liam
Liam 2022-11-16
@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
Matt J 2022-11-17
I can't provide my raw data because it has identifiable donor health information.
But why 2 .mat files instead of 1.

请先登录,再进行评论。

采纳的回答

Matt J
Matt J 2022-11-16
编辑:Matt J 2022-11-16
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
  3 个评论
Matt J
Matt J 2022-11-18
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
Matt J 2022-11-20
编辑:Matt J 2022-11-20
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

请先登录,再进行评论。

更多回答(1 个)

Alex Sha
Alex Sha 2022-11-17
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]
  6 个评论
Matt J
Matt J 2022-11-20
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
Alex Sha 2022-11-20
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.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Dimensionality Reduction and Feature Extraction 的更多信息

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by