How to use Nlinfit for a function with two independent variables?

Hi here is my data and code and I am trying to predict parameters for a function with two independent variables but Nlinfit is giving me error.

clc

clear %all; % Clear the workspace.

close all; % Close all figures.

format compact

%% Read in data

data =readmatrix('Flexible_BUC.xlsx');

x1=[7.50000000000000

7.50000000000000

7.70000000000000

5

5

5

5

5

5

5

5

5

7.50000000000000

7.50000000000000

7.50000000000000

7.50000000000000

7.50000000000000

6

6

6

6

6

6

3.75000000000000

3.75000000000000

3.75000000000000

3.75000000000000

3.75000000000000

8

8

8

8

8

8

7.50000000000000

7.50000000000000

7.50000000000000

8

8

8];

x2=[0.00153000000000000

0.00522000000000000

0.000189000000000000

3.73000000000000e-06

3.73000000000000e-06

1.17000000000000e-05

8.66000000000000e-05

1.17000000000000e-05

1.51000000000000e-05

2.99000000000000e-05

7.00000000000000e-05

7.92000000000000e-05

6.06000000000000e-05

0.000163000000000000

0.000656000000000000

0.000818000000000000

0.00129000000000000

9.78000000000000e-06

0.000161000000000000

0.000183000000000000

0.000204000000000000

0.000297000000000000

0.000343000000000000

0.00705000000000000

0.00850000000000000

0.0233000000000000

0.0250000000000000

0.0267000000000000

0.00125000000000000

0.00245000000000000

0.00391000000000000

0.00878000000000000

0.00111000000000000

0.00122000000000000

3.23000000000000e-05

4.26000000000000e-05

4.53000000000000e-05

4.13000000000000e-05

8.24000000000000e-05

8.33000000000000e-05];

yobs=[5.95833333300000

0.300000000000000

0.0625000000000000

0.111111111000000

0.213809289000000

0.140625000000000

0.651315789000000

0.351694915000000

12.0555555600000

0.626846311000000

0.555555556000000

0.136363636000000

6.38793103400000

0.233051458000000

0.540000000000000

0.0277777780000000

1.05555555600000

0.113636364000000

0.0933323590000000

0.352272727000000

3.20833333300000

0.897435897000000

1.17046404700000

1.41666666700000

1.79545454500000

1.15384615400000

1.85576923100000

10.8333333300000

0.848684211000000

6.18835443000000

0.767441860000000

0.527777778000000

1.54872306000000

0.337691494000000

1.08333333300000

1.87477002000000

1.81654734900000

1.97222222200000

0.550000000000000

0.340020401000000];

xm=[x1 x2];

%% Initial parameter guesses

C1=0.25;

C2=0.73;

beta0(1)=C1; %initial guess beta 1

beta0(2)=C2; %initial guess beta 2

p=length(beta0); %p = # parameters

%% define function to be used for inverse problem

fINV=@Project_funcINV;

%fnameINV=@forderexpINV;

[beta,resids,J,COVB,mse] = nlinfit(xm,yobs,fINV,beta0);

beta

%% Functions

function y=Project_funcINV(beta0,x1,x2)

c2s=@(x)-2.40874-39.748*(1+x).^-2.856;

y=100./(1+exp(-beta0(1).*c2s(x1)+(beta0(2).*c2s(x1).*log10(100.*x2))));

end

Someone can please help, I would appreciate it.

Answer 1

[Star Strider]

You are using the correct approach with:

xm=[x1 x2];

In the function, refer to ‘x1’ as ‘xm(:,1)’ and ‘x2’ as ‘xm(:,2)’ , passing ‘xm’ as the independent variable to ‘Project_funcINV’. I made those changes, and added a fitnlm call to display the statistics, and provided a plot of the data and the fit to it (as a line plot).

Try this —

x1=[7.50000000000000
7.50000000000000
7.70000000000000
5
5
5
5
5
5
5
5
5
7.50000000000000
7.50000000000000
7.50000000000000
7.50000000000000
7.50000000000000
6
6
6
6
6
6
3.75000000000000
3.75000000000000
3.75000000000000
3.75000000000000
3.75000000000000
8
8
8
8
8
8
7.50000000000000
7.50000000000000
7.50000000000000
8
8
8];
x2=[0.00153000000000000
0.00522000000000000
0.000189000000000000
3.73000000000000e-06
3.73000000000000e-06
1.17000000000000e-05
8.66000000000000e-05
1.17000000000000e-05
1.51000000000000e-05
2.99000000000000e-05
7.00000000000000e-05
7.92000000000000e-05
6.06000000000000e-05
0.000163000000000000
0.000656000000000000
0.000818000000000000
0.00129000000000000
9.78000000000000e-06
0.000161000000000000
0.000183000000000000
0.000204000000000000
0.000297000000000000
0.000343000000000000
0.00705000000000000
0.00850000000000000
0.0233000000000000
0.0250000000000000
0.0267000000000000
0.00125000000000000
0.00245000000000000
0.00391000000000000
0.00878000000000000
0.00111000000000000
0.00122000000000000
3.23000000000000e-05
4.26000000000000e-05
4.53000000000000e-05
4.13000000000000e-05
8.24000000000000e-05
8.33000000000000e-05];
yobs=[5.95833333300000
0.300000000000000
0.0625000000000000
0.111111111000000
0.213809289000000
0.140625000000000
0.651315789000000
0.351694915000000
12.0555555600000
0.626846311000000
0.555555556000000
0.136363636000000
6.38793103400000
0.233051458000000
0.540000000000000
0.0277777780000000
1.05555555600000
0.113636364000000
0.0933323590000000
0.352272727000000
3.20833333300000
0.897435897000000
1.17046404700000
1.41666666700000
1.79545454500000
1.15384615400000
1.85576923100000
10.8333333300000
0.848684211000000
6.18835443000000
0.767441860000000
0.527777778000000
1.54872306000000
0.337691494000000
1.08333333300000
1.87477002000000
1.81654734900000
1.97222222200000
0.550000000000000
0.340020401000000];
xm=[x1 x2];
%% Initial parameter guesses
C1=0.25;
C2=0.73;
beta0(1)=C1; %initial guess beta 1
beta0(2)=C2; %initial guess beta 2
p=length(beta0); %p = # parameters
%% define function to be used for inverse problem
fINV=@Project_funcINV;
%fnameINV=@forderexpINV;
[beta,resids,J,COVB,mse] = nlinfit(xm,yobs,fINV,beta0);
beta

beta = 1×2

1.3585 0.1749

mdl = fitnlm(xm,yobs,fINV,beta0)                            % ADDED

mdl = 
Nonlinear regression model:
    y ~ Project_funcINV(b,X)

Estimated Coefficients:
          Estimate      SE       tStat       pValue  
          ________    _______    ______    __________

    b1     1.3585     0.13278    10.231    1.8021e-12
    b2    0.17495     0.11064    1.5812       0.12212


Number of observations: 40, Error degrees of freedom: 38
Root Mean Squared Error: 2.81
R-Squared: -0.0164,  Adjusted R-Squared -0.0431
F-statistic vs. zero model: 7.46, p-value = 0.00185

figure                                                      % ADDED
stem3(x1, x2, yobs, 'filled')
hold on
plot3(x1, x2, Project_funcINV(beta,xm), '-r')
hold off

%% Functions
function y=Project_funcINV(beta0,xm)
c2s=@(x)-2.40874-39.748*(1+x).^-2.856;
y=100./(1+exp(-beta0(1).*c2s(xm(:,1))+(beta0(2).*c2s(xm(:,1)).*log10(100.*xm(:,2)))));
end

The fit is reasonably good, although ‘beta(2)’ may not be significnatly different from zero.

.

[Faizan Lali]

Thank you very much. It worked.

[Star Strider]

As always, my pleasure!

[Faizan Lali]

Moreover, if we ignore the likely outliers, what do you say about the function and parameter estimates? Because I am getting the zero in confidence intervals for Beta 1, is there any way to estimate significantly ?

[Star Strider]

I am not surprised about the confidence interval being essentially 0, considering that the p value that fitnlm calculates for it is , likely the significance value you want. Statistics calculations for it are likely close to the eps value as the result.

[Faizan Lali]

Gotcha.... Thank you

[Star Strider]

As always, my pleasure!

[Faizan Lali]

Hi,

I have a quick question extension to previous one. Once I calculate the 95% confidence and prediction bands for depndent variable it gives 40x1 column vector for each. So I am curious how to plot confidence interval and y_pred ? and agaisnt which "x" I should plot them ?

Is there any way out to plot them in 3 D ?

here is my code I am using to calculate confidence and prediction bands.

You help will be beneficial. Thank you!

ci=nlparci(beta, resids,J)

%R is the correlation matrix for the parameters, sigma is the standard error vector

[R,sigma]=corrcov(COVB);

R

R_2=R.*R;

%Rsq1 = 1 - sum((yobs - ypred).^2)/sum((yobs - mean(yobs)).^2);

sigma

relerr=sigma./beta' %relative error (coefficient of variance) for each parameter

%% Confidence and prediction intervals for the dependent variable

[ypred, delta] = nlpredci(fnameINV,xm,beta,resids,J,0.05,'on','curve'); %confidence band for regression line

[ypred, deltaob] =nlpredci(fnameINV,xm,beta,resids,J,0.05,'on','observation');%prediction band for individual points

yspan=range(ypred)% total span of ypred

relrmse=rmse/yspan % ratio of rmse vs. yspan

%simultaneous confidence bands for regression line

CBu=ypred+delta;

CBl=ypred-delta;

[Star Strider]

You can use any ‘x’ values you want, providing that they span the same bounds as the original ‘x’ values. This means that you can use the linspace function to create them (although I do not recommend that with this data set). If you use fitnlm, you can use the appropriate predict function.

x1=[7.50000000000000
7.50000000000000
7.70000000000000
5
5
5
5
5
5
5
5
5
7.50000000000000
7.50000000000000
7.50000000000000
7.50000000000000
7.50000000000000
6
6
6
6
6
6
3.75000000000000
3.75000000000000
3.75000000000000
3.75000000000000
3.75000000000000
8
8
8
8
8
8
7.50000000000000
7.50000000000000
7.50000000000000
8
8
8];
x2=[0.00153000000000000
0.00522000000000000
0.000189000000000000
3.73000000000000e-06
3.73000000000000e-06
1.17000000000000e-05
8.66000000000000e-05
1.17000000000000e-05
1.51000000000000e-05
2.99000000000000e-05
7.00000000000000e-05
7.92000000000000e-05
6.06000000000000e-05
0.000163000000000000
0.000656000000000000
0.000818000000000000
0.00129000000000000
9.78000000000000e-06
0.000161000000000000
0.000183000000000000
0.000204000000000000
0.000297000000000000
0.000343000000000000
0.00705000000000000
0.00850000000000000
0.0233000000000000
0.0250000000000000
0.0267000000000000
0.00125000000000000
0.00245000000000000
0.00391000000000000
0.00878000000000000
0.00111000000000000
0.00122000000000000
3.23000000000000e-05
4.26000000000000e-05
4.53000000000000e-05
4.13000000000000e-05
8.24000000000000e-05
8.33000000000000e-05];
yobs=[5.95833333300000
0.300000000000000
0.0625000000000000
0.111111111000000
0.213809289000000
0.140625000000000
0.651315789000000
0.351694915000000
12.0555555600000
0.626846311000000
0.555555556000000
0.136363636000000
6.38793103400000
0.233051458000000
0.540000000000000
0.0277777780000000
1.05555555600000
0.113636364000000
0.0933323590000000
0.352272727000000
3.20833333300000
0.897435897000000
1.17046404700000
1.41666666700000
1.79545454500000
1.15384615400000
1.85576923100000
10.8333333300000
0.848684211000000
6.18835443000000
0.767441860000000
0.527777778000000
1.54872306000000
0.337691494000000
1.08333333300000
1.87477002000000
1.81654734900000
1.97222222200000
0.550000000000000
0.340020401000000];
xm=[x1 x2];
%% Initial parameter guesses
C1=0.25;
C2=0.73;
beta0(1)=C1; %initial guess beta 1
beta0(2)=C2; %initial guess beta 2
p=length(beta0); %p = # parameters
%% define function to be used for inverse problem
fINV=@Project_funcINV;
%fnameINV=@forderexpINV;
[beta,resids,J,COVB,mse] = nlinfit(xm,yobs,fINV,beta0);
beta

beta = 1×2

1.3585 0.1749

mdl = fitnlm(xm,yobs,fINV,beta0)                            % ADDED

mdl = 
Nonlinear regression model:
    y ~ Project_funcINV(b,X)

Estimated Coefficients:
          Estimate      SE       tStat       pValue  
          ________    _______    ______    __________

    b1     1.3585     0.13278    10.231    1.8021e-12
    b2    0.17495     0.11064    1.5812       0.12212


Number of observations: 40, Error degrees of freedom: 38
Root Mean Squared Error: 2.81
R-Squared: -0.0164,  Adjusted R-Squared -0.0431
F-statistic vs. zero model: 7.46, p-value = 0.00185

[ypredc, deltac] = nlpredci(@Project_funcINV,xm,beta,resids,'Jacobian',J,'Alpha',0.05,'SimOpt','on','PredOpt','curve'); %confidence band for regression line
% [ypredob, deltaob] =nlpredci(fnameINV,xm,beta,resids,J,0.05,'on','observation');%prediction band for individual points

CBu=ypredc+deltac;
CBl=ypredc-deltac;

figure                                                      % ADDED
stem3(x1, x2, yobs, 'filled')
hold on
plot3(x1, x2, Project_funcINV(beta,xm), '-r')
plot3(x1, x2, [CBu CBl], '--r')
hold off

%% Functions
function y=Project_funcINV(beta0,xm)
c2s=@(x)-2.40874-39.748*(1+x).^-2.856;
y=100./(1+exp(-beta0(1).*c2s(xm(:,1))+(beta0(2).*c2s(xm(:,1)).*log10(100.*xm(:,2)))));
end

.

[Faizan Lali]

Thank you very much. It helped alot.

[Star Strider]

As always, my pleasure!

[Alex Sha]

In theory, the following result is a little better (the result provided by Star Strider is a local solution):

Sum Squared Error (SSE): 297.811503388026
Root of Mean Square Error (RMSE): 2.72860542854782
Correlation Coef. (R): 0.527048050632212
R-Square: 0.277779647675215

Parameter	Best Estimate   
---------	-------------   
b1   	9.64813154107489
b2   	20.8981735230572

[Faizan Lali]

How you parameter estimates are so high?

Which code did you use?

[Walter Roberson]

Alex uses a commercial program named 1stOpt that does some very nice optimization. Sometimes I am able to improve a little on his results, but not usually, and when I do manage then it is only after a couple of days of continuous computations.