n=2;% размерность случайного вектора
N=1000;%количество наблюдений (объем выборки)
U = lognrnd(0,1,n,N);%генерация Х по закону распределения Релея
cov_u=eye(n,n); %ковариационная матрица
Dx=1;
Dy=1;
cov_ksi=[Dx 0;0 Dy]; %ввод ковариационной матрицы
M_ksi=[0;0] ;

A=zeros(n,n);
%проверка размерности cov_ksi и M_ksi
if (size(cov_ksi) ~= n) | (size(M_ksi) ~= n)
 error('Размерность матрицы сov_ksi или M_ksi не совпадает с n');
end
for i=1:n
 if det(cov_ksi(1:i,1:i)) <= 0
 error('Матрица сov_ksi не положительно определенна');
 end
end

%вычисление элементов матрицы A
for i=1:n
 for j=1:i
 sum=0;
 for k=1:(j-1)
 sum=sum+A(i,k)*A(j,k)*cov_u(k,k);
 end
 if i==j
 A(i,j)=sqrt((cov_ksi(i,j)-sum)/cov_u(j,j));
 else
 A(i,j)=(cov_ksi(i,j)-sum)/(A(j,j)*cov_u(j,j));
 end
 end
end

%построение случайного вектора ksi
for i=1:N
 ksi(:,i)=A*U(:,i)+M_ksi;
end
%транспонирование матрицы случайных векторов ksi
ksi_t=ksi';
disp('Матрица преобразований A');
A;
disp('Выборка входных векторов U');
U;
disp('Выходные векторы ksi');
ksi;
disp('Исходн. матрица cov_ksi');
cov_ksi;
disp('Матрица cov_ksi (проверка правильности преобразований)');
A*cov_u*A';
disp('Выборочная ков.матрица');
test_cov_ksi=ksi'/(n-1);

%построение гистограмм для входных и выходных векторов
figure(1);
hist(U',20);
xlabel('Интервалы');
ylabel(sprintf('Количество из %g',N));
title('Гистограмма каждой из компоненты вход. двум.СВ(x1,x2)');

figure(2);
hist(ksi',50);
xlabel('Интервалы');
ylabel(sprintf('Количество из %g',N));
title('Гистограмма каждой из компоненты выход. двум.СВ(x1,x2)')

figure(3);
subplot(2,1,1);
plot(U');
xlabel('Номер эксперемента,n');
ylabel('Пара (x1,x2) n-го экспер');
title('Входной двум.СВ(x1,x2)')
figure(3);
subplot(2,1,2);
plot(ksi');
xlabel('Номер эксперемента,n');
ylabel('Пара (x1,x2) n-го экспер');
title('Выходной двум.СВ(x1,x2)')

% %плотность вероятности двумерной величины
m=500;% размер массива
x_ksi=sqrt(Dx);
y_ksi=sqrt(Dy);
r=(mean(ksi(1)*ksi(2))-mean(ksi(1))*mean(ksi(2)))/(x_ksi*y_ksi);%коэфициент корреляции
mx=M_ksi(1);
my=M_ksi(2);
z=zeros(m);
test_z=zeros(m);
for i=1:m
 x(i)=i/30-9;
 for j=1:m
 y(j)=j/30-9;

 ch=1/(2*pi*x_ksi*y_ksi*sqrt(1-r^2));
 power_1=-1/(2*(1-r^2));
 power_2=((x(i)-mx)^2)/(x_ksi)^2;
 power_3=(2*r*(x(i)-mx)*(y(j)-my))/(x_ksi*y_ksi);
 power_4=((y(j)-my)^2)/(y_ksi)^2;
 rez_exp=exp(power_1*(power_2-power_3+power_4));
 z(i,j)= rez_exp*ch;
 end
end

x_ksi=test_cov_ksi(1,1);
y_ksi=test_cov_ksi(2,2);
r=(mean(ksi(1)*ksi(2))-mean(ksi(1))*mean(ksi(2)))/(x_ksi*y_ksi);%коэфициент корреляции
mx=M_ksi(1);
my=M_ksi(2);

for i=1:m
 x(i)=i/30-9;
 for j=1:m
 y(j)=j/30-9;

 ch=1/(2*pi*x_ksi*y_ksi*sqrt(1-r^2));
 power_1=-1/(2*(1-r^2));
 power_2=((x(i)-mx)^2)/(x_ksi)^2;
 power_3=(2*r*(x(i)-mx)*(y(j)-my))/(x_ksi*y_ksi);
 power_4=((y(j)-my)^2)/(y_ksi)^2;
 rez_exp=exp(power_1*(power_2-power_3+power_4));
 test_z(i,j)= rez_exp*ch;
 end
end

figure();
mesh(x,y,z);
xlabel('x1');
ylabel('x2');
title('p(x1,x2)')
title('Плотность (теор) вероятности двум.СВ (x1,x2)');

figure();
mesh(x,y,test_z);
xlabel('x1');
ylabel('x2');
title('p(x1,x2)')
title('Плотность (экспр) вероятности двум.СВ (x1,x2)');

figure();
[C,h]=contourf(x,y,z);
clabel(C,h);
xlabel('x1');
ylabel('x2');
title('Контурный график (теория)');

figure();
[C,h]=contourf(x,y,test_z);
clabel(C,h);
xlabel('x1');
ylabel('x2');
title('Контурный график (эксперимент)');