How to obtain the standard deviation of the fitting parameters?

26 次查看(过去 30 天)
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)));

采纳的回答

Star Strider
Star Strider 2024-3-21
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)));
.
  6 个评论
Qili Hu
Qili Hu 2024-3-22
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
Star Strider 2024-3-22
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.
.

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Calculus 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by