from sympy.abc import y,x
from sympy import integrate, simplify

fxy = x + y                 # joint density
fy = integrate(fxy,(x,0,1)) # marginal density
fx = integrate(fxy,(y,0,1)) # marginal density

h = (3*y+2)/(6*y+3) # conditional expectation
LHS=integrate((x - h)**2 *fxy, (x,0,1),(y,0,1)) # from the definition
RHS=integrate( (x)**2 *fx, (x,0,1)) - integrate( (h)**2 *fy, (y,0,1)) # using Pythagorean theorem
print simplify(LHS-RHS)==0