//main
clear all;
close all;
clc

A = load('/Users/marialukasova/Downloads/CD_ИАД_2019-2/ЛР_2/Data_Lab2/data4.txt');

x = A(:,1);
y = A(:,2);
sigma = A(:,3);
figure;
plot(x,y,'black'); %A(:,1)-- удаление первого столбца A(:,2) -- второго
title ('экспирементальные данные')
xlabel('Argument x')
ylabel('Function y')

dop = A(:,1);
dop = dop(:,1);
s = size(dop);
a = [2, 3];

%relej
[teta1, r_rel, j_rel, cov_rel] = nlinfit(x, y, @rel, a);
y_rel = rel(teta1, x);


%vejbull
[teta2, r_vej, j_vej, cov_vej] = nlinfit(x, y, @vej, a);
y_vej = vej(teta2, x);

%gamma
[teta3, r_gam, j_gam, cov_gam] = nlinfit(x, y, @gam, a);
figure;
plot (x, y, 'black');
title ('Аппроксимация гамма');
xlabel ('X');
ylabel ('Y');
hold on
y_gam = gam(teta3, x);
plot (x, y_gam, 'g', x, y_rel, 'r', x, y_vej, 'b');
hold off

ksi2_rel = ksi2(y_rel, y, sigma, 200-2-1, s);
ksi2_vej = ksi2(y_vej, y, sigma, 200-2-1, s);
ksi2_gam = ksi2(y_gam, y, sigma, 200-2-1, s);

res_rel = res(y_rel, y, sigma);
res_vej = res(y_vej, y, sigma);
res_gam = res(y_gam, y, sigma);

figure;
subplot(3, 1, 1);
plot(x,res_rel, x, zeros(1,200),'r');
title('График взвешенных остатков Релей');
xlabel('X');
ylabel('R');
subplot(3, 1, 2);
plot(x,res_vej, x, zeros(1,200),'r');
title('График взвешенных остатков Вейбулл');
xlabel('X');
ylabel('R');
subplot(3, 1, 3);
plot(x,res_gam, x, zeros(1,200),'r');
title('График взвешенных остатков гамма');
xlabel('X');
ylabel('R');

N=1:100;
a_rel = avtokor(res_rel);
a_vej = avtokor(res_vej);
a_gam = avtokor(res_gam);
figure;
subplot(3, 1, 2);
plot(N, a_vej,N,zeros(1,100),'r');
title ('Автокорреляционная функция для Вейбула')
xlabel ('k')
ylabel ('A(k)')
subplot(3, 1, 1);
plot(N, a_rel,N,zeros(1,100),'r');
title ('Автокорреляционная функция для Релея')
xlabel ('k')
ylabel ('A(k)')
subplot(3, 1, 3);
plot(N, a_gam,N,zeros(1,100),'r');
title ('Автокорреляционная функция для гаммы')
xlabel ('k')
ylabel ('A(k)')


dov_rel1 = nlparci(teta1, r_rel, 'jacobian',j_rel, 'alpha', 0.68);
dov_vej1 = nlparci(teta2, r_vej, 'jacobian',j_vej, 'alpha', 0.68);
dov_gam1 = nlparci(teta3, r_gam, 'jacobian',j_gam, 'alpha', 0.68);

diag_cov_rel = diag(cov_rel);
dov_rel2 = [teta1(1) - tinv(0.68, 199)*sqrt(diag_cov_rel(1)) teta1(1) + tinv(0.68, 199)*sqrt(diag_cov_rel(1)); teta1(2) - tinv(0.68, 199)*sqrt(diag_cov_rel(2)) teta1(2) + tinv(0.68, 199)*sqrt(diag_cov_rel(2))];


diag_cov_vej = diag(cov_vej);
dov_vej2 = [teta2(1) - tinv(0.68, 199)*sqrt(diag_cov_vej(1)) teta2(1) + tinv(0.68, 199)*sqrt(diag_cov_vej(1)); teta2(2) - tinv(0.68, 199)*sqrt(diag_cov_vej(2)) teta2(2) + tinv(0.68, 199)*sqrt(diag_cov_vej(2))];

diag_cov_gam = diag(cov_gam);
dov_gam2 = [teta3(1) - tinv(0.68, 199)*sqrt(diag_cov_gam(1)) teta3(1) + tinv(0.68, 199)*sqrt(diag_cov_gam(1)); teta3(2) - tinv(0.68, 199)*sqrt(diag_cov_gam(2)) teta3(2) + tinv(0.68, 199)*sqrt(diag_cov_gam(2))];


//relej
function y=rel(sig,x)
y=x./(sig(1)^2).*exp(-(x.^2)/2/(sig(2)^2));



//vejbull
function y=vej(beta, x)
y = beta(1).*beta(2).*(x.^(beta(2)-1)).*exp(-beta(1).*(x.^beta(2)));



//gamma
function result = gam(a, x)
result = x.^a(1).*exp(-x/a(2))/(a(2).^(a(1)+1)*gamma(a(1)+1));



//ksi_kv
function ksi2 = ksi2(f, y, sigma,v, s)
ks=0;
for k=1:s
    ks = ks + ((f(k)-y(k))/sigma(k)).^2;
end
ksi2 = ks/v;


//ostatki
function r = res(f, y, sigma)

for k=1:200
    r(k) = (f(k) - y(k))/sigma(k);
end




//avtokor
function a = avtokor(R)
dop = 0;
for i = 1:200
    dop = dop + R(i).^2;
end

for k = 1:100
    dop2 = 0;
    for i = 1:(200 - k+1)
        dop2= dop2 + R(i) .* R(i+k-1);
    end
    a(k)=(dop2./dop).*200./(200 - k + 1);
end