How to obtain the standard deviation of the fitting parameters?
15 次查看(过去 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)));
0 个评论
采纳的回答
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)
CovMtx = CovB
ci = nlparci(beta1,Rv,'Covar',CovB)
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 个评论
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'
ci=nlparci(beta1,Rv,'Covar',CovB) % PArameter Confidence Intervals
tci=tinv([0.025 0.975],12-2); % Inverse 't' Distribution
sigma=((ci-betav)./tci)*sqrt(12) % Calculate Standard Deviations
fprintf('Beta(1) = %8.4f\t\tStandatrd Deviation = %8.4f\nBeta(2) = %8.4f\t\tStandatrd Deviation = %8.4f\n',[betav sigma(:,1)].')
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 Center 和 File Exchange 中查找有关 Calculus 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!