Independence of $(X,Y)$ with joint PDF $f_{X,Y}(x,y)=8xy$ on $ 0I know that as a general property, if
$f_{X,Y}(x,y)=g(x)h(y)$, then we say X and Y are independent random variables.
However, I am having trouble to accept this statement. Take the case:
$f_{X,Y}(x,y)=8xy, 0<y<x<1$
Then, $f_{X}(x)=4x$ and $f_{Y}(y)=4y$, where $f_{X}(x).f_{Y}(y)\ne f_{X,Y}(x,y)$ 
Maybe I am getting something wrong over the limits of integration, I often trip on simple things... Or maybe, in this case, since the variables are dependent in the limits of integration, then they cannot be independent. But I havent seen any book point that as a requirement... So I suppose I got something wrong in my calculations. 
If someone could please clarify this for me, I would be grateful.
 A: 
Take the case: $f_{X,Y}(x,y)=8xy, 0<y<x<1$. 

Yet again an excellent example of the fact that densities should include restrictions on their domain.
Recall that any density of the couple of random variables $(X,Y)$ should be a function $f:\mathbb R\times\mathbb R\to\mathbb R_+$ such that what-you-know holds. In the present case, $f:\mathbb R\times\mathbb R\to\mathbb R_+$, $(x,y)\mapsto8xy$ is (obviously) not a density. Rather, a density of $(X,Y)$ is  $f:\mathbb R\times\mathbb R\to\mathbb R_+$, $(x,y)\mapsto8xy\mathbf 1_{0<y<x<1}$ (for example $f(1,2)=0$, not $16$).
Now the result you mention is correct:

The random variables $X$ and $Y$ with density $f$ are independent if and only if there exist $g$ and $h$ such that $f(x,y)=g(x)h(y)$ for (almost) every $(x,y)$ in $\mathbb R\times\mathbb R$.

Which does not hold for the density $f$ in the example.
This remark is also useful when computing marginals. In general,
$$
f_X(x)=\int_\mathbb Rf(x,y)\mathrm dy,
$$
hence in the present case,
$$
f_X(x)=\mathbf 1_{0<x<1}\int_\mathbb R8xy\mathbf 1_{0<y<x}\mathrm dy=x\mathbf 1_{0<x<1}\int_0^x8y\mathrm dy=4x^3\mathbf 1_{0<x<1}.
$$
Note the absence of cases and the automaticity of the computations, thanks to adequate notations.
To sum up:

When a joint distribution is given by its PDF, a détour by the joint CDF is useless (and frankly often cumbersome) provided one uses the true PDF, which should include indicator functions if need be.

A: The complete, formal, heavy, official, no-joking, definition of independence is: Two random variables $X$ and $Y$ are independent if every function of $X$ is independent of every function of $Y$... Thankfully, it suffices to state what @Zelareth stated, namely that they are independent iff their joint probability density function can be written as the product of their marginal densities (what if they don't have densities?)
$$ $$
As for your example, it is one of the tricky ones. Indeed, the variables are not independent. To derive their density functions, the safest approach (but visibly longer) is to go all the way up to the joint distribution function, then back down to the marginal distribution functions and then to the marginal density functions. But this is safe, and illuminating, especially when the support sets of the variables are "linking" the variables together. In general, when the supports are not plus-minus infinity, care must always be exercised, even if they appear disjoint.
So by definition the joint distribution function is $F(x,y)= P(\min x\le X\le x, \min y\le Y\le y)$. In our case $\min X = y$, while $\min Y=0$. Using dummy variables of integration ($u$ for $x$ and $w$ for $y$), we have 
$$F(x,y)=\int_{w=0}^y \int_{u=w}^xf_{XY}(u,w)dudw=8\int_{w=0}^yw \int_{u=w}^xududw=8\int_{w=0}^yw \frac12\left(x^2-w^2\right)dw=$$
$$4x^2\int_{w=0}^yw dw -4\int_{w=0}^yw^3dw=2x^2y^2-y^4 \;,\;  0<y<x<1\qquad [1]$$
One critical detail to remember in the above, is that any lower limits of integration that involve $x$ or $y$ must be expressed in terms of the dummy variables of integration, while the upper integration limits must be kept written in terms of the actual variables $X$ and $Y$.
We turn now to the marginal distribution functions. By definition
$$F(x) = \lim_{y\to \max y}F(x,y)$$
In our case (this is another critical detail) $\max y = x$. And this holds although the marginal density of $y$ will have support $(0,1)$. Intuitively, we haven't yet "separated" the variables, so we must still respect their interrelation.  Substituting in $F(x,y)$ we obtain
$$F(x) = \lim_{y\to x}\left(2x^2y^2-y^4\right) = 2x^4-x^4 = x^4 \;,\; x\in (0,1)\qquad [2]$$
Now that we have ousted $Y$, the variable $X$ can behave as though $Y$ doesn't exist, and so its support is $(0,1)$. The marginal density of $X$ is the derivative:
$$f_X(x)=\frac{d}{dx}F(x) = 4x^3\;,\; x\in (0,1) \qquad [3]$$
You can verify that it integrates to unity over its support.
For the $Y$ variable we have analogously
$$F(y) = \lim_{x\to \max x}F(x,y)$$
In our case  $\max x = 1$. Substituting in $F(x,y)$ we obtain
$$F(y) = \lim_{x\to 1}\left(2x^2y^2-y^4\right) = 2y^2-y^4 \;,\; y\in (0,1)\qquad [4]$$
and the density is of $Y$ is:
$$f_Y(y)=\frac{d}{dy}F(y) = 4y-4y^3\;,\; y \in (0,1) \qquad [5]$$
It too integrates to unity. As you can see, the product of the marginal densities has nothing to do with the joint density, so the variables are not independent. I would suggest you work the conditional distributions and densities, to complete the example.
A: I think you're slightly confused with this characterization of independence. If we assume that $X$ and $Y$ have density functions $f_X$ and $f_Y$ and joint distribution $f_{X,Y}$, then $X$ and $Y$ are independent iff $f_{X,Y} = f_x \cdot f_y$. This is not the same as the property you listed above, the $g(x)$ and $h(y)$ must be the marginal density functions and not just any function of $x$ and $y$.
