mashlukashova

Untitled

Mar 18th, 2021
1,024
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. //main
  2. clear all;
  3. close all;
  4. clc
  5.  
  6. A = load('/Users/marialukasova/Downloads/CD_ИАД_2019-2/ЛР_2/Data_Lab2/data4.txt');
  7.  
  8. x = A(:,1);
  9. y = A(:,2);
  10. sigma = A(:,3);
  11. figure;
  12. plot(x,y,'black'); %A(:,1)-- удаление первого столбца A(:,2) -- второго
  13. title ('экспирементальные данные')
  14. xlabel('Argument x')
  15. ylabel('Function y')
  16.  
  17. dop = A(:,1);
  18. dop = dop(:,1);
  19. s = size(dop);
  20. a = [2, 3];
  21.  
  22. %relej
  23. [teta1, r_rel, j_rel, cov_rel] = nlinfit(x, y, @rel, a);
  24. y_rel = rel(teta1, x);
  25.  
  26.  
  27. %vejbull
  28. [teta2, r_vej, j_vej, cov_vej] = nlinfit(x, y, @vej, a);
  29. y_vej = vej(teta2, x);
  30.  
  31. %gamma
  32. [teta3, r_gam, j_gam, cov_gam] = nlinfit(x, y, @gam, a);
  33. figure;
  34. plot (x, y, 'black');
  35. title ('Аппроксимация гамма');
  36. xlabel ('X');
  37. ylabel ('Y');
  38. hold on
  39. y_gam = gam(teta3, x);
  40. plot (x, y_gam, 'g', x, y_rel, 'r', x, y_vej, 'b');
  41. hold off
  42.  
  43. ksi2_rel = ksi2(y_rel, y, sigma, 200-2-1, s);
  44. ksi2_vej = ksi2(y_vej, y, sigma, 200-2-1, s);
  45. ksi2_gam = ksi2(y_gam, y, sigma, 200-2-1, s);
  46.  
  47. res_rel = res(y_rel, y, sigma);
  48. res_vej = res(y_vej, y, sigma);
  49. res_gam = res(y_gam, y, sigma);
  50.  
  51. figure;
  52. subplot(3, 1, 1);
  53. plot(x,res_rel, x, zeros(1,200),'r');
  54. title('График взвешенных остатков Релей');
  55. xlabel('X');
  56. ylabel('R');
  57. subplot(3, 1, 2);
  58. plot(x,res_vej, x, zeros(1,200),'r');
  59. title('График взвешенных остатков Вейбулл');
  60. xlabel('X');
  61. ylabel('R');
  62. subplot(3, 1, 3);
  63. plot(x,res_gam, x, zeros(1,200),'r');
  64. title('График взвешенных остатков гамма');
  65. xlabel('X');
  66. ylabel('R');
  67.  
  68. N=1:100;
  69. a_rel = avtokor(res_rel);
  70. a_vej = avtokor(res_vej);
  71. a_gam = avtokor(res_gam);
  72. figure;
  73. subplot(3, 1, 2);
  74. plot(N, a_vej,N,zeros(1,100),'r');
  75. title ('Автокорреляционная функция для Вейбула')
  76. xlabel ('k')
  77. ylabel ('A(k)')
  78. subplot(3, 1, 1);
  79. plot(N, a_rel,N,zeros(1,100),'r');
  80. title ('Автокорреляционная функция для Релея')
  81. xlabel ('k')
  82. ylabel ('A(k)')
  83. subplot(3, 1, 3);
  84. plot(N, a_gam,N,zeros(1,100),'r');
  85. title ('Автокорреляционная функция для гаммы')
  86. xlabel ('k')
  87. ylabel ('A(k)')
  88.  
  89.  
  90. dov_rel1 = nlparci(teta1, r_rel, 'jacobian',j_rel, 'alpha', 0.68);
  91. dov_vej1 = nlparci(teta2, r_vej, 'jacobian',j_vej, 'alpha', 0.68);
  92. dov_gam1 = nlparci(teta3, r_gam, 'jacobian',j_gam, 'alpha', 0.68);
  93.  
  94. diag_cov_rel = diag(cov_rel);
  95. 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))];
  96.  
  97.  
  98. diag_cov_vej = diag(cov_vej);
  99. 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))];
  100.  
  101. diag_cov_gam = diag(cov_gam);
  102. 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))];
  103.  
  104.  
  105. //relej
  106. function y=rel(sig,x)
  107. y=x./(sig(1)^2).*exp(-(x.^2)/2/(sig(2)^2));
  108.  
  109.  
  110.  
  111. //vejbull
  112. function y=vej(beta, x)
  113. y = beta(1).*beta(2).*(x.^(beta(2)-1)).*exp(-beta(1).*(x.^beta(2)));
  114.  
  115.  
  116.  
  117. //gamma
  118. function result = gam(a, x)
  119. result = x.^a(1).*exp(-x/a(2))/(a(2).^(a(1)+1)*gamma(a(1)+1));
  120.  
  121.  
  122.  
  123. //ksi_kv
  124. function ksi2 = ksi2(f, y, sigma,v, s)
  125. ks=0;
  126. for k=1:s
  127.     ks = ks + ((f(k)-y(k))/sigma(k)).^2;
  128. end
  129. ksi2 = ks/v;
  130.  
  131.  
  132. //ostatki
  133. function r = res(f, y, sigma)
  134.  
  135. for k=1:200
  136.     r(k) = (f(k) - y(k))/sigma(k);
  137. end
  138.  
  139.  
  140.  
  141.  
  142. //avtokor
  143. function a = avtokor(R)
  144. dop = 0;
  145. for i = 1:200
  146.     dop = dop + R(i).^2;
  147. end
  148.  
  149. for k = 1:100
  150.     dop2 = 0;
  151.     for i = 1:(200 - k+1)
  152.         dop2= dop2 + R(i) .* R(i+k-1);
  153.     end
  154.     a(k)=(dop2./dop).*200./(200 - k + 1);
  155. end
  156.  
RAW Paste Data