Текст программы:
clear all;
nx=41;
ny=10;
f_read=fopen('data.txt', 'r');
[x,nxx]=fscanf(f_read,['%f' ' '],[nx,1]);
[y,nyy]=fscanf(f_read,['%f' '\n'],[ny,inf]);
f=fopen('results.txt','w');
x=x';
y=sort(y);
y_sr=mean(y);
y_disp=var(y);
q=0.95;
k=3.379; % толерантный множитель для n=10, P=0.95, Q=0.95
t_lim=[y_sr-k.*sqrt(y_disp);y_sr+k.*sqrt(y_disp)];
dov_int=[y_sr-sqrt(y_disp./ny)*tinv((1+q)/2,(ny-1)); y_sr-sqrt(y_disp./ny)*tinv((1+q)/2,(ny-1))];
krit_znach=0.1133; % для n=10 k=41
krit_kochr=max(y_disp)/sum(y_disp);
fprintf(f,'Статистика критерия Кочрена: %f\n', krit_kochr);
fprintf(f,'Критическое значение: %f\n', krit_znach);
p=-1;
q=0;
fprintf(f, '\nq F критическое F\n');
while(p<0)
X=ones(nx,q+1);
for i=0:q
X(:,i+1)=(x').^i;
end
a=(X'*X)^(-1)*X'*y_sr';
R2=ny*(sum(y_disp/nx))^(-1)*(X*a-y_sr')'*(X*a-y_sr');
F=R2/(nx-q-1);
F_krit=finv(0.95,nx-q-1,nx-1);
fprintf(f, '%d %f %f\n',q,F,F_krit);
if(F<=F_krit)
p=q;
else
q=q+1;
end
end
fprintf(f, '\nКовариационная матрица при q: \n');
for j=1:nx
s2(j)=1/(ny-1)*sum((y(:,j)-y_sr(j)).^2);
end
Kov=sum(s2)/nx/ny*(X'*X)^(-1);
for i=1:size(Kov,1)
for j=1:size(Kov,2)
fprintf(f, ['%0.3e' ' '],Kov(i,j));
end
fprintf(f, '\n');
end
cond_kov=cond(Kov);
fprintf(f, 'Cond= %d \n',cond_kov);
Kor=zeros(q+1,q+1);
fprintf(f, '\nКорреляционная матрица при q:\n');
for i=1:q+1
for j=1:q+1
Kor(i,j)=Kov(i,j)/sqrt(Kov(i,i)*Kov(j,j));
fprintf(f, ['%0.3e' ' '],Kor(i,j));
end
fprintf(f, '\n');
end
fprintf(f, '\nq F критическое F\n');
q1=0;
while(q1<nx)
X1=ones(nx,q1+1);
for i=0:q1
X1(:,i+1)=(x').^i;
end
a2=(X1'*X1)^(-1)*X1'*y_sr';
R2_2=ny*(sum(y_disp/nx))^(-1)*(X1*a2-y_sr')'*(X1*a2-y_sr');
F=R2_2/(nx-q1-1);
F_krit=finv(0.95,nx-q1-1,nx-1);
fprintf(f, '%d %f %f\n',q1,F,F_krit);
q1=q1+1;
end
Kov2=sum(s2)/nx/ny*(X1'*X1)^(-1);
cond_kov2=cond(Kov2);
fprintf(f, '\nпри k-1: Cond= %d \n',cond_kov2);
Kor2=zeros(nx,nx);
for i=1:k
for j=1:k
Kor2(i,j)=Kov2(i,j)/sqrt(Kov2(i,i)*Kov2(j,j));
end
end
for i=1:q+1
polynom_q(i)=a((q+1)+1-i);
end
for i=1:nx
polynom_k(i)=a2(nx+1-i);
end
figure(1);
plot(x,y_sr,'ko',x,t_lim(1,:),'k--',x,t_lim(2,:),'k--',x,dov_int(1,:),'g--',x,dov_int(2,:),'g--',...
x,polyval(polynom_q,x),'r*-',x,polyval(polynom_k,x),'b*-');
legend('Математическое ожидание','толерантные пределы','толерантные пределы','границы доверительных интервалов',...
'границы доверительных интервалов','полином степени q','полином степени k-1',0);
figure(2);
regr=regress(y_sr',[ones(nx,1) x']);
robust=robustfit(x,y_sr');
polynom=polyfit(x,y_sr,1);
ridg=ridge(y_sr',[ones(nx,1) x'],0);
plot(x,polyval(regr(end:-1:1),x),'r-', x,polyval(robust(end:-1:1),x),'g-',x,polyval(polynom,x),'b-', x,polyval(ridg(end:-1:1),x),...
'c-',x,y_sr,'m*',x,y,'k.');
legend('метод R-Square','робастная регрессия',...
'полиномиальная регрессия','ридж-регрессия','математическое ожидание','исходные точки');
title('Аппроксимация прямой линией');
%polytool(x,y_sr)
figure(3);
plot(x,polyval(polyfit(x,y_sr,12),x),'r-',x,y_sr,'bo',x,y,'k.');
title('Полиномиальная аппроксимация');
legend('аппроксимирующая линия','математическое ожидание','исходные точки');
x1=min(x):((x(2)-x(1))/10):max(x);
figure(4);
inter1=interp1(x,y_sr,x);
inter2=interp1(x,y_sr,x1,'cubic');
erm=pchip(x,y_sr,x1);
spl=spline(x,y_sr,x1);
plot(x,inter1,'g-', x1,inter2,'r-',x1,erm,'b-', x1,spl,'c-',x,y_sr,'m*',x,y,'k.');
title('Кусочная полиномиальная аппроксимация');
legend('линейная','кубическая','полиномы Эрмита','сплайн','математическое ожидание','исходные точки');
figure(5);
b0=ones(1,7);
b=nlinfit(x',y_sr',@func,b0);
plot(x,func(b,x),x,y_sr,'ro',x,y,'k.')
title('Нелинейная аппроксимация');
legend('аппроксимирующая линия','математическое ожидание','исходные точки');
fprintf(f,'\n Коэффициенты для нелинейной аппроксимации:\n');
fprintf(f,'%f\n',b);