% Valores de una v.a. X que toma 4 valores entre 100 y 200 (vector columna)
lx = 4;
X = sort(randi([100,200],lx,1)); % imponemos orden ascendiente
% Valores de una v.a. Y que toma 5 valores entre 50 y 150 (vector fila)
ly = 5;
Y = sort(randi([50,150],1,ly));
% Definir la matriz de la distribución conjunta del vector (X,Y)
P = randi(30,lx,ly);
P = (1/sum(sum(P))) * P ;
% Distribuciones marginales
PY = sum(P); % suma por columnas
PX = sum(transpose(P));
% 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,length(XC)) * P * reshape(YC,length(YC),1) ;
% 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 = reshape(X,lx,1) * ones(1,ly) + ones(lx,1) * reshape(Y,1,ly);
% Una pequen~a comprobacion
if any( size(SS) ~= size(P) )
   disp('SS y P deberian tener las mismas dimensiones');
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-10;
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 = '); disp(P);
   fprintf('Var(X) = '); disp(VX);
   fprintf('Var(Y) = '); disp(VY);
   fprintf('Cov(X,Y) = '); disp(CovXY);
   fprintf('Var(X+Y) = '); disp(VS);
   fprintf('Var(X) + Var(Y) + 2 * Cov(X,Y) = '); disp(VX+VY+2*CovXY);
end