% Valores de una v.a. X que toma valores entre 100 y 2000 (vector fila)
clear all;
X = sort(unique(randi([100,2000],1,8))); % imponemos orden ascendiente estricto
lx = length(X);
% Valores de una v.a. Y que toma valores entre 50 y 1500 (vector fila)
Y = sort(unique(randi([50,1500],1,10)));
ly = length(Y);
% Definir la matriz de la distribucion conjunta del vector (X,Y)
P = randi([100,300],lx,ly); % frecuencias absolutas
P = (1/sum(sum(P))) * P ;   % frecuencias relativas (probabilidades)
% Distribuciones marginales
PX = sum(P,2); % para cada fila, suma por columnas
PY = sum(P,1); % para cada columna, suma por filas
% Esperanzas
EX = dot( PX, X);
EY = dot( PY, Y);
% Variables aleatorias centradas
XC = X - EX * ones(size(X));
YC = Y - EY * ones(size(Y));
% Varianzas
VX = dot( PX, XC.^2 );
VY = dot( PY, YC.^2 );
% Covarianza del vector (X,Y)
% (aseguramos compatibilidad de dimensiones mediante 'reshape')
CovXY = reshape(XC,1,lx) * P * reshape(YC,ly,1) ;
CorrXY = CovXY / sqrt(VX * VY);
% Ahora calcularemos la suma, su distribucion y su varianza directamente
% En primer lugar, creamos una matriz con las dimensiones de P que
%  contiene las respectivas sumas
SS = zeros(lx,ly);
% En notacion matricial:
% SS = reshape(X,lx,1) * ones(1,ly) + ones(lx,1) * reshape(Y,1,ly);
% Alternativamente, lo hacemos elemento a elemento:
for i = 1 : lx
    for j = 1 : ly
        SS(i,j) = X(i) + Y(j);
    end
end
% Ahora pasamos las sumas y su distribucion a vectores fila
S = reshape(SS,1,prod(size(SS)));
PS = reshape(P,1,prod(size(P)));
% Esperanza de la suma
ES = dot(PS, S);
% Suma centrada
SC = S - ES * ones(size(S));
% Varianza de la suma
VS = dot(PS, SC.^2 );
% Se debe cumplir que VS = VX + VY + 2 * CovXY, vamos a comprobarlo
tol = 1e-6;
if abs( VX + VY + 2 * CovXY - VS ) > tol
   disp('Ha habido algun error en los calculos de VS');
else
   disp('El calculo de la varianza de la suma coincide con los dos metodos');
   fprintf('X = '); disp(X');
   fprintf('Y = '); disp(Y);
   fprintf('P = \n'); disp(P);
   fprintf('Var(X) = '); disp(VX);
   fprintf('Var(Y) = '); disp(VY);
   fprintf('Cov(X,Y) = '); disp(CovXY);
   fprintf('Corr(X,Y) = '); disp(CorrXY);
   fprintf('Var(X+Y) = '); disp(VS);
   fprintf('Var(X) + Var(Y) + 2 * Cov(X,Y) = '); disp(VX+VY+2*CovXY);
end