function c = corr_from_joint_distribution(P,x_values, y_values)

assert(ndim(P)==2);
assert(length(x_values)==size(P,1));
assert(length(y_values)==size(P,2));

P = P./sum(P(:));
Px = sum(P,2);
Py = sum(P,1);

mX = dot(x_values, Px);
mY = dot(y_values, Py);
XY = reshape(x_values,length(x_values),1)*reshape(y_values,1,length(y_values));
mXY = dot(XY(:),P(:));

varX = dot(x_values.^2, Px) - mX^2;
varY = dot(y_values.^2, Py) - mY^2;

covXY = mXY - mX*mY;
c = covXY/sqrt(varX*varY);

end