%% 
P = [0.02 0.10 0.08; 0.08 0.32 0.40];
x = [1 3]; y = [2 4 5];

pX = sum(P,2)'
pY = sum(P)

EX = x*pX'

VarX = ((x-EX).^2)*pX'

EY = y*pY'

VarY = ((y-EY).^2)*pY'

[X Y] = meshgrid(x,y); X = X'; Y = Y';

PXYl6 = sum(sum([X.*Y<6].*P))
PXeY = sum(sum([X==Y].*P))

EXY = sum(sum(X.*Y.*P))

EX3Y2 = sum(sum((X-3).*(Y-2).*P))

EXY311Y238Y = sum(sum( (X).*(Y.^3 - 11*(Y.^2) + 38*Y).*P ))

CovXY = EXY - EX*EY

rhoXY = CovXY/(sqrt(VarX*VarY))