Azzera filtri
Azzera filtri

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

18 visualizzazioni (ultimi 30 giorni)
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 Commenti
Liam
Liam il 16 Nov 2022
@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 il 17 Nov 2022
I can't provide my raw data because it has identifiable donor health information.
But why 2 .mat files instead of 1.

Accedi per commentare.

Risposta accettata

Matt J
Matt J il 16 Nov 2022
Modificato: Matt J il 16 Nov 2022
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 Commenti
Matt J
Matt J il 18 Nov 2022
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 il 20 Nov 2022
Modificato: Matt J il 20 Nov 2022
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

Accedi per commentare.

Più risposte (1)

Alex Sha
Alex Sha il 17 Nov 2022
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 Commenti
Matt J
Matt J il 20 Nov 2022
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 il 20 Nov 2022
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.

Accedi per commentare.

Prodotti


Release

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by