How to obtain the standard deviation of the fitting parameters?
24 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Qili Hu
el 21 de Mzo. de 2024
Comentada: Star Strider
el 22 de Mzo. de 2024
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 comentarios
Respuesta aceptada
Star Strider
el 21 de Mzo. de 2024
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 comentarios
Star Strider
el 22 de Mzo. de 2024
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.
.
Más respuestas (0)
Ver también
Categorías
Más información sobre Calculus en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!