Гамма - распределение и критерий согласия Пирсона

ymax2=max(wy2);
ymin2=min(wy2);
 
k=round(n^0.5); % число интервалов для построения гистограммы
d=ymax2-ymin2/k; % ширина каждого интервала
 
clear Tabl % очистили таблицу результатов
Tabl(:,1)=[1:k]'; % номера интервалов
 
% Запись значений с шагом d
yml2(1)=0;
for i=2:1:k
    yml2(i)=yml2(i-1)+d;
end;
yml2=yml2';
Tabl(:,2)=yml2;
 
ymp2(1)=d;
 
for i=2:1:k
ymp2(i)=ymp2(i-1)+d;
end;
 
ymp2=ymp2';
 
Tabl(:,3)=ymp2;
 
 
nj=hist(wy2,k);
Tabl(:,4)=nj';
 
bor=[Tabl(:,2);Tabl(end,3)]; % все границы интервалов
 
mx=mean(wy2);
sx=std(wy2);
lam=abs(1/mx); % параметр показательного распределения
a=mx-sx*3^0.5; 
b=mx+sx*3^0.5; % параметры равномерного распределения
sig=abs(mx)*(2/pi)^0.5; % параметр Рэлеевского распределения
 
pro=[normcdf(bor,mx,sx),expcdf(bor,1/lam),unifcdf(bor,a,b),raylcdf(bor,sig)];
 
nd=4;
Tabl(:,5)=pro(2:end,nd)-pro(1:end-1,nd); % вероятности попадания pj
 
Tabl(:,6)=n*Tabl(:,5); % теоретическое число попаданий npj
Tabl(:,7)=(Tabl(:,4)-Tabl(:,6)).^2./Tabl(:,6); % относит.квадр.отклонение
 
disp('Сводная таблица результатов')
disp(' j aj bj nj pj npj (nj-npj)^2/npj')
fprintf('%2.0f%12.5f%12.5f%6.0f%12.5f%12.5f%12.5f\n',Tabl')
 
ResTabl=Tabl(1,1:6); % взяли первую строку
 
for k1=2:k, % берём остальные строки таблицы
 
 if ResTabl(end,6)<5, % предыдущее npj<5 => будем суммировать
 
 ResTabl(end,3)=Tabl(k1,3); % новая правая граница интервала
 
 ResTabl(end,4:6)=ResTabl(end,4:6)+Tabl(k1,4:6); % суммируем nj, pj, npj
 
 else % предыдущее npj>=5 => будем дописывать строку
 
 ResTabl=[ResTabl;Tabl(k1,1:6)]; % дописываем строку
 
 end
 
end
 
if ResTabl(end,6)<5, % последнее npj<5
 
 ResTabl(end-1,3)=ResTabl(end,3); % новая правая граница интервала
 
 ResTabl(end-1,4:6)=ResTabl(end-1,4:6)+ResTabl(end,4:6); % nj, pj, npj
 
ResTabl=ResTabl(1:end-1,:); % отбросили последнюю строку
 
end
 
kn=size(ResTabl,1); % число объединённых интервалов
 
ResTabl(:,1)=[1:kn]'; % новые номера интервалов
 
ResTabl(:,7)=(ResTabl(:,4)-ResTabl(:,6)).^2./ResTabl(:,6); % (nj-npj)^2/npj
 
disp('Сгруппированная сводная таблица результатов')
 
disp(' j aj bj nj pj npj (nj-npj)^2/npj')
 
fprintf('%2.0f%12.5f%12.5f%6.0f%12.5f%12.5f%12.5f\n',ResTabl')
 
hi2=sum(ResTabl(:,7)); % сумма элементов последнего столбца - статистика Пирсона
 
fprintf('Статистика Пирсона chi2=%10.5f\n',hi2)
 
pz=0.1; % выбрали уровень значимости
 
fprintf('Задаём уровень значимости p=%5.4f\n',pz)
 
chi2pz=chi2inv(1-pz,kn-3); % квантиль распределения Пирсона
 
fprintf('Квантиль chi2-распределения Пирсона chi2(%5.2f,%d)=%10.5f\n',1-pz,kn-3, chi2pz)
 
if hi2<=chi2pz,
 
 disp('Распределение подобрано верно, так как chi2<=chi2(1-p,k-3)')
 
else
 
 disp('Распределение подобрано неверно, так как chi2>chi2(1-p,k-3)')
 
end

d=(ymax2-ymin2)/k; % ширина каждого интервала