funkce=A*Sin[x*y]+A;
intervalX={x,0,1.};
intervalY={y,0,1.};
(*kovariance[funkce_]:=Module[{res,margY,margX,stredniX,stredniY,stredniXY,kovariance},*)
res=Solve[Integrate[Integrate[funkce,intervalX],intervalY]==1,A][[1]]
margY=Integrate[funkce/.res,intervalX];
margX=Integrate[funkce/.res,intervalY];

stredniX=Integrate[margX*x,intervalX];
stredniY=Integrate[margY*y,intervalY];
 stredniXY=Integrate[Integrate[x*y*funkce/.res,intervalX],intervalY ];
kovariance=stredniX*stredniY-stredniXY