from sympy import Matrix, Symbol, exp, pi, simplify, integrate from sympy import stats, sqrt, oo, use from sympy.abc import y,x sigma_x = Symbol('sigma_x',positive=True) sigma_y = Symbol('sigma_y',positive=True) sigma_xy = Symbol('sigma_xy',real=True) fyy = stats.density(stats.Normal('y',0,sigma_y))(y) R = Matrix([[sigma_x**2, sigma_xy], [sigma_xy,sigma_y**2]]) fxy = 1/(2*pi)/sqrt(R.det()) * exp((-Matrix([[x,y]])*R.inv()* Matrix([[x],[y]]))[0,0]/2 ) fcond = simplify(fxy/fyy) u=Symbol('u',positive=True) # define positive variable fcond2=fcond.subs(sigma_x**2*sigma_y**2-sigma_xy**2,u) # substitute as hint to integrate g=simplify(integrate(fcond2*x,(x,-oo,oo))) # evaluate integral gg=g.subs( u,sigma_x**2 *sigma_y**2 - sigma_xy**2 ) # substitute back in use( gg, simplify,level=2) # simplify exponent term