How to obtain the standard deviation of the fitting parameters?

I have a question about how to obtain the standard deviation of the fitting parameters according to the following program codes.

x=[0 5 10 15 20 30 45 60 75 90 105 120];

y=[0 3.87 4.62 4.98 5.21 5.40 5.45 5.50 5.51 5.52 5.54 5.53];

plot(x,y,'bo');

hold on

pause(0.1);

beta0=[39,0.002];

% syms n t

% fun=@(beta,t) beta(1)*(1-6/(pi^2)*symsum((1./n.^2).*exp(-beta(2)*(n.^2).*t),n,1,Inf));

% betafit = nlinfit(x,y,fun,beta0);

beta1=beta0;

delta = 1e-8; % desired objective accuracy

R0=Inf; % initial objective function

for K=1:10000

fun=@(beta,t) beta(1)*(1-6/(pi^2)*sum((1./(1:K)'.^2).*exp(-beta(2)*((1:K)'.^2).*t),1));

[betafit,R] = nlinfit(x,y,fun,beta1);

R = sum(R.^2);

if abs(R0-R)<delta

break;

end

beta1=betafit;

R0 = R;

end

plot(x,fun(betafit,x),'.-r');

xlabel('x');

ylabel('y');

legend('experiment','model');

title(strcat('\beta=[',num2str(betafit),'];----stopped at--','K=',num2str(K)));

Answer 1

[Star Strider]

The most meaningful measure of the parameters is the confidence interval. You could possibly recover some information from the covariance matrix ‘CovB’, however that would take some effort.

x=[0 5 10 15 20 30 45 60 75 90 105 120];
y=[0 3.87 4.62 4.98 5.21 5.40 5.45 5.50 5.51 5.52 5.54 5.53];
plot(x,y,'bo');
hold on
pause(0.1);
beta0=[39,0.002];
% syms n t
% fun=@(beta,t) beta(1)*(1-6/(pi^2)*symsum((1./n.^2).*exp(-beta(2)*(n.^2).*t),n,1,Inf));
% betafit = nlinfit(x,y,fun,beta0);
beta1=beta0;
delta = 1e-8; % desired objective accuracy
R0=Inf; % initial objective function
% K=1:10000;
for K=1:10000
    fun=@(beta,t) beta(1)*(1-6/(pi^2)*sum((1./(1:K)'.^2).*exp(-beta(2)*((1:K)'.^2).*t),1));
    [betafit,Rv,J,CovB,MSE] = nlinfit(x,y,fun,beta1);
    R = sum(Rv.^2);
    if abs(R0-R)<delta
        break;
    end
    beta1=betafit;
    R0 = R;
end
fprintf('\nbeta = %8.4f\nbeta = %8.4f\n\n',beta1)

beta =   5.4900
beta =   0.1374

CovMtx = CovB

CovMtx = 2×2

1.0e-03 * 0.4686 -0.0534 -0.0534 0.0206

ci = nlparci(beta1,Rv,'Covar',CovB)

ci = 2×2

5.4418 5.5382 0.1273 0.1475

plot(x,fun(betafit,x),'.-r');
xlabel('x');
ylabel('y');
legend('experiment','model', 'Location','best');
title(strcat('\beta=[',num2str(betafit),'];----stopped at--','K=',num2str(K)));

.

[Qili Hu]

Dear Strider

I appreciate for your warm work earnestly. I do not quite understand the codes you write. I want to konw the standard deviation of two paramenters qe and D in the following equation. Can you help me to write it in detail? Thank you very much.

the standard deviation of the fitting parameters

[Star Strider]

This calculation is not usually performed, since the standard deviation of a parameter in a multi-parametric system is not important.

The standard error of the estimate is with σ denoting the standard deviation, and the confidence intervals multiply it by the critical values of the t distribution for the problem being considered, . The confidence intervals use the inverse t distribution with n-2 degrees -of-freedom (since this estimates 2 parameters), so one way of calculating it could be —

tci = tinv([0.025 0.975], 12-2)

tci = 1×2

-2.2281 2.2281

beta =   [5.4900
          0.1374];                              % Parameter Estimates
ci = [5.4418    5.5382
      0.1273    0.1475];                        % Confidence Intervals

sigma = ((ci-beta)./tci) * sqrt(12)             % Estimated Standard Deviations

sigma = 2×2

0.0749 0.0749 0.0157 0.0157

This is the best that I can do to estimate the standard deviations, based on assumptions regarding confidence intervals using the standard errors of the mean for vectors. I cannot claim that it is accurate, since I have never done this calculation before, nor have I ever seen it done.

.

[Qili Hu]

Dear Strider

I appreciate for solvingthis question for me earnestly. Thank you very much again.

Bist wishes.

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.

[Qili Hu]

Dear Strider

According to the codes you write, I add them to my codes. But, it do not work during running. When I change sigma=((ci-beta)./tci)*sqrt(12) to sigma=((ci-beta1)./tci)*sqrt(12), the running result is not consistent with you. It seems that beta do not appear in the codes. So, I hope you can help me to check the following codes. Thank you very much.

x=[0 5 10 15 20 30 45 60 75 90 105 120];

y=[0 3.87 4.62 4.98 5.21 5.40 5.45 5.50 5.51 5.52 5.54 5.50];

plot(x,y,'ro','MarkerFaceColor','r');

hold on

beta0=[5,0.1];

beta1=beta0;

delta=1e-8;

R0=Inf;

for n=1:15000

fun=@(beta,t) beta(1)*(1-6/(pi^2)*sum((1./(1:n)'.^2).*exp(-beta(2)*((1:n)'.^2).*t),1));

[betafit,Rv,J,CovB,MSE]=nlinfit(x,y,fun,beta1);

R=sum(Rv.^2);

if abs(R0-R)<delta

break;

end

beta1=betafit;

R0=R;

end

CovMtx=CovB;

ci=nlparci(beta1,Rv,'Covar',CovB);

tci=tinv([0.025 0.975],12-2);

sigma=((ci-beta)./tci)*sqrt(12);

plot(x,fun(betafit,x),'-r');

xlabel('t (min)');

ylabel('qt (mg/g)');

title(strcat('\beta=[',num2str(betafit),'];----stopped at--','n=',num2str(n)));

[Star Strider]

It seems that you only copied part of my code.

When I provided the rest of it (slightly modified from the earlier version), it works as expected —

x=[0 5 10 15 20 30 45 60 75 90 105 120];
y=[0 3.87 4.62 4.98 5.21 5.40 5.45 5.50 5.51 5.52 5.54 5.50];
plot(x,y,'ro','MarkerFaceColor','r');
hold on
beta0=[5,0.1];
beta1=beta0;
delta=1e-8;
R0=Inf;
for n=1:15000
   fun=@(beta,t) beta(1)*(1-6/(pi^2)*sum((1./(1:n)'.^2).*exp(-beta(2)*((1:n)'.^2).*t),1));
   [betafit,Rv,J,CovB,MSE]=nlinfit(x,y,fun,beta1);
   R=sum(Rv.^2);
   if abs(R0-R)<delta
      break;
   end
   beta1=betafit;
   R0=R;
end
% CovMtx=CovB;
betav = beta1(:)                                % Convert 'beta1' To The column Vector 'betav'

betav = 2x1

5.4856 0.1379

ci=nlparci(beta1,Rv,'Covar',CovB)               % PArameter  Confidence Intervals

ci = 2x2

5.4387 5.5326 0.1280 0.1478

tci=tinv([0.025 0.975],12-2);                   % Inverse 't'  Distribution
sigma=((ci-betav)./tci)*sqrt(12)                % Calculate Standard Deviations

sigma = 2x2

0.0731 0.0731 0.0154 0.0154

fprintf('Beta(1) = %8.4f\t\tStandatrd Deviation = %8.4f\nBeta(2) = %8.4f\t\tStandatrd Deviation = %8.4f\n',[betav sigma(:,1)].')

Beta(1) =   5.4856		Standatrd Deviation =   0.0731
Beta(2) =   0.1379		Standatrd Deviation =   0.0154

plot(x,fun(betafit,x),'-r');
xlabel('t (min)');
ylabel('qt (mg/g)');
title(strcat('\beta=[',num2str(betafit),'];----stopped at--','n=',num2str(n)));

I added the fprintf call to be certain that there are no ambiguities.

.