Confidence and prediction bands using nlsqcurvefit

I am trying to fit the model first by nlinfit and then lsqcurvefit but in both cases I am getting the non significant parameter estimation for one of the parameter as its confidence intervals conatains zero. There are two limitations of my model;

1. I can not change the model or can not transform any of the data.
2. Parameters can not be less than zero.

Working on the 2nd point I used lsqcurvefit with bound limits for parameters as

lb = [0.001,0.1];
ub = [1,1.5];

but unfortunately the 2nd parameter does not change its value from lower bound. In other words while estimating by lsqcurvefit the function is always giving the value of beta2 equal to corresponding lower bound value.

Also I want to plot Confidence & prediction bands for the estimates. I did it with nlinfit but I am not sure how to do it with lsqcurvefit.

Here is my code I Don't know where I am making mistake. Data file is attached.

data =readmatrix('BUC_Matlab_Input.xlsx');
x1=data(:,2);
x2=data(:,3);
yobs=data(:,4);
xm=[x1 x2];
%% Initial parameter guesses
% C1=0.32;
% C2=0.27;
C1=1.31;
C2=2.1585;
beta0(1)=C1; %initial guess yo
beta0(2)=C2; %initial guess kr
p=length(beta0); %p = # parameters
%% define the function that will be used for the forward problem 
fnameFOR = @Project_func;
ypred_b=fnameFOR(beta0,x1,x2);
%% lsqcurvefit returns parameters, residuals, Jacobian (sensitivity %coefficient matrix),
%how to find covariance matrix, and mean square error ??
fnameINV=@Project_funcINV;
lb = [0.001,0.1];
ub = [1,1.5];
[beta,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(fnameINV,beta0,xm,yobs,lb,ub);
beta
% Confidence intervals using nlparci
conf = nlparci(beta,residual,'jacobian',jacobian);
% Confidence interval calculated manually
lastwarn(''); %Clear warning memory
alpha = 0.05; %significance level
df = length(residual) - numel(beta); %degrees of freedom
crit = tinv(1-alpha/2,df);       %critical value
covm = inv(jacobian'*jacobian) * var(residual); %covariance matrix
[~, warnId] = lastwarn; %detect if inv threw 'nearly singular' warning.
covmIdx = sub2ind(size(covm),1:size(covm,1),1:size(covm,2));  %indices of the diag of covm
CI = nan(numel(beta),2);
if ~strcmp(warnId, 'MATLAB:nearlySingularMatrix')
    CI(:,1) = beta - crit * sqrt(covm(covmIdx));
    CI(:,2) = beta + crit * sqrt(covm(covmIdx));
end
%% Confidence and Prediction bands
% How to calculate CBs and PBs ?? I have included the code below which I am using
% in nlinfit estimation.
%% Confidence and prediction intervals for the dependent variable using nlinfit (This part is just for reference)
%nonlinear regression confidence intervals-- 'on' means simultaneous
%bounds; 'off' is for nonsimultaneous bounds; must use 'curve' for
%regression line, 'observation' for prediction interval
[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;
%simultaneous prediction bands for regression line
PBu=ypred+deltaob;
PBl=ypred-deltaob;
%% Functions
% Forward Function
function y=Project_func(beta0,x1,x2)
%tm=[x1 x2];
c2s=@(x)-2.40874-(39.748*(1+x).^-2.856);
y=100./(1+exp((-beta0(1).*2.*c2s(x1))+(beta0(2).*c2s(x1).*log10(100.*x2))));
end
% Inverse Function
function y=Project_funcINV(beta0,xm)
c2s=@(x)-2.40874-39.748*(1+x).^-2.856;
y=100./(1+exp(-beta0(1).*2.*c2s(xm(:,1))+(beta0(2).*c2s(xm(:,1)).*log10(100.*xm(:,2)))));
end

Answer 1

[Matt J]

Covariance methods are not valid when you have parameter bound constraints. Since it's not a time-consuming fit, I would just use some sort of Monte Carlo resampling to determine the variability of the beta estimates:

data =readmatrix('BUC_Matlab_Input.xlsx');
x1=data(:,2);
x2=data(:,3);
yobs=data(:,4);
xm=[x1 x2];
%% Initial parameter guesses
% C1=0.32;
% C2=0.27;
C1=1.31;
C2=2.1585;
beta0(1)=C1; %initial guess yo
beta0(2)=C2; %initial guess kr
p=length(beta0); %p = # parameters

%% lsqcurvefit returns parameters, residuals, Jacobian (sensitivity %coefficient matrix),
%how to find covariance matrix, and mean square error ??
fnameINV=@Project_funcINV;
lb = [0.001,0.1];
ub = [1,1.5];
[beta,resnorm,residual,exitflag] = lsqcurvefit(fnameINV,beta0,xm,yobs,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.

sig=std(residual);
ypred=fnameINV(beta,xm);

BETAS=monteCarlo(fnameINV,beta,xm,ypred,lb,ub,sig); %Monte Carlo samples


std(BETAS)

ans = 1×2

0.0619 0.0820

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

function BETAS=monteCarlo(fnameINV,beta,xm,ypred,lb,ub,sig)

  N=100; %number of random samples

  BETAS=nan(N,numel(beta));
  opts=optimoptions('lsqcurvefit','Display','none');

  for i=1:N
   yobs=ypred+sig*randn(size(ypred));
   BETAS(i,:) = lsqcurvefit(fnameINV,beta,xm,yobs,lb,ub,opts);
  end

end

[Faizan Lali]

Thank you @Matt J, I appreciate it, but I have another code, where I am using ga to minimise std error and it is giving me beta values as [0.76185 0.01624], which are generally acceptable for this model. If you want I can attach that codes, but its kind a mess.

Moreover, my main focus it to determine and plot the Confidence and prediction bands for these estimates.

[Torsten]

Read again Matt's answer:

Covariance methods are not valid when you have parameter bound constraints.

[Matt J]

Thank you @Matt J, I appreciate it

You're quite welcome, but if this answers your question, please Accept-click the answer.

[Faizan Lali]

Once I use these betas which you suggested getting from montecarlo, my functions prediction is not quite right... Is there other way to optimize this funcion?

[Matt J]

The optimization is done in the first call to lsqcurvefit. The BETAS are just re-estimates to allow you to assess the statistical uncertainty of the estimates. They are deliberately wrong.

[Faizan Lali]

Thank you so much @Matt J, I got it.

[Faizan Lali]

@Matt J But in first call the estimation of the betas remains same, as beta 2 is only picking ub the lb value 0.1 in this case. it will do the same once we give such kind constraints to betas estimation? I am just kinda curious now.

[Matt J]

I don't know what you mean by "the betas remain the same". If you look at the 100 different re-estimations in BETAS, you will see that neither beta(1) nor beta(2) remain the same.

[Faizan Lali]

@Matt J I was talking about this;

[beta,resnorm,residual,exitflag] = lsqcurvefit(fnameINV,beta0,xm,yobs,lb,ub);

but I got it now, that how monteCarlo is working.

Can we find rmse or relerror in this setup?

[Matt J]

Do you know ground truth?

[Faizan Lali]

No.. What do you mean by ground truth?

[Matt J]

Calculating rmse and relerror require that you know the true value of beta.