Transforming any probability density over $R^N$ to any other probability density If $p_1(x)$ and $p_2(x)$ are two arbitrary probability density functions defined over $R^N$, i.e. both are non-negative and properly normalized so that
$$ \int d^N x \, p_1(x) = 1 = \int d^N x \, p_2(x) \,, $$ 
then does there exist a change of variable $y=y(x)$ that maps one density to the other? In other words, does there exist a $y(x)$ such that
$$ d^N x \, p_1(x) = d^N x \Big|\frac{\partial y}{\partial x} \Big|\, p_2(y(x)) \,,$$
where $\Big|\frac{\partial y}{\partial x} \Big|$ is the determinant of the Jacobian due to the change of variables?
 A: There exists a differentiable change of variable when both densities are everywhere nonzero.
To see why, let's reduce the question to a simpler one through a series of steps.  The following process could be expressed as a (rather complicated) explicit formula in terms of $p_1$ and $p_2,$ but the formula would be opaque: it's good to know the demonstration is constructive and can even be implemented in software.
First, to streamline the exposition let's take $N=2:$ the generalization to higher dimensions will be clear.  We are thereby considering random variables $(X_1,Y_1)$ with density $p_1$ and $(X_2,Y_2)$ with density $p_2.$
Second, let the marginal (cumulative) distribution functions of $X_1$ and $Y_1$ be $F_1$ and $G_1,$ respectively. The map
$$\phi_1:(X_1,Y_1) \to (F_1(X_1),G_1(Y_1)) = (U_1,V_1)$$
(sometimes called the Probability Integral Transformation) defines a distribution $C_1$ on the unit square $[0,1]^N=[0,1]\times[0,1].$ $C_1$ is called a copula: see Nelsen.
Explicitly, for any $(x,y)\in \mathbb{R}^2,$
$$C_1(F_1(x), G_1(y)) = \int_{-\infty}^x \int_{-\infty}^y p_1(x,y)\,\mathrm{d}y\,\mathrm{d}x.$$
Analogously define $F_2,$ $G_2,$ $\phi_2,$ $U_2,$ $V_2,$ and $C_2$ for density $p_2.$  If we can transform $C_1$ to $C_2$ via a diffeomorphism $\psi$ of $[0,1]^2,$ then we can transform $p_1$ to $p_2$ by applying the composition $\phi_2^{-1}\circ \psi \circ \phi_1.$
Third, without any loss of generality we may take $C_2$ to be the independence copula $\Pi$ with constant density $1$ on $[0,1]^2:$ $$\Pi(u,v)=uv\text{ for }(u,v)\in[0,1]^2.$$  If we can transform any density to $\Pi$ invertibly, then we need only transform $C_1$ to $\Pi$ and follow that with the inversion of the transformation from $C_2$ to $\Pi.$
Finally, let $C$ be either $C_1$ or $C_2:$ that is, it's an arbitrary copula for a continuous joint distribution with nonzero density $p.$  The map
$$\psi: (u,v) \to \left(u, \frac{C(u,v)}{u}\right)$$
transforms $C$ to $\Pi$ because

*

*$$0 \le \frac{C(u,v)}{u} \le \frac{C(u,1)}{u} = \frac{u}{u} = 1$$ shows $\psi: [0,1]^2\to [0,1]^2,$


*$\psi$ is invertible because $p\ne 0$ implies that for any $(x,y)\in[0,1]^2$ there is a unique $(u,v)$ for which $\psi(u,v)=(x,y),$ and


*$$\Pr(U\le u,\ V\le v) = C(u,v) = u\, \frac{C(u,v)}{u} = \Pi\left(u, \frac{C(u,v)}{u} \right)$$ shows $\psi$ converts $C$ into $\Pi,$ QED.
At step (2) we needed at least one of the densities to be everywhere nonzero.  This suggests there may be counterexamples to the general proposition.  One counterexample is afforded by any density $p$ on $[0,1]^2$ that equals zero on $[1/2,1]\times[0,1/2]$ (such as the uniform density on the complement of that square).  Any invertible map $\psi$ must produce a density that is zero in some neighborhood of $(3/4,1/4)$ and therefore cannot be the uniform density on $[0,1]^2.$

If you had in mind the more flexible situation of an arbitrary change of variable (although preferably almost everywhere differentiable), then we could eliminate the restriction to nonzero densities.
As before, reduce the problem to transforming one copula into another, where each copula is determined by a density.  The density defines an (infinite) quadtree.  The root node is the rectangle $[0,1]^2.$  Split that rectangle along the largest coordinate $x_1$ for which $C(x_1,1)=1/2:$ the left rectangle is $[0,x_1]\times[0,1]$ and the right rectangle is its complement $(x_1,1]\times[0,1].$  Then, separately split the two resulting rectangles along the largest coordinates $y_1$ for which $C(x_1,y_1)=1/4$ (for the left rectangle) and $y_2$ for which $C(1,y_2)=3/4$ (for the right rectangle).  This partitions $[0,1]^2$ into four disjoint rectangles, each of probability $1/4:$ they form the children of the root node.  Adopt a convention to put them in a sequence (such as lower left, upper left, lower right, and upper right).
Proceed recursively to partition each child (in order) into four disjoint rectangles, each with one-quarter the probability of its parent.  At each stage of this recursion all the terminal nodes are given a definite sequence.  See my answer at https://gis.stackexchange.com/a/31879/664 for an illustration of this process with large discrete datasets. Here is an illustration showing the rectangles produced after approximately five iterations, along with points representing the density (think of them as a stratified sample of that density):

Two such quadtrees determine a mapping $\psi$ from $[0,1]^2$ to itself. Approximate $\psi$ by considering the two sequences of $4^n$ rectangles established after level $n$ of the recursion.  Let $R_i$ be one of those rectangles for the first quadtree and $S_i$ the corresponding rectangle for the second quadtree.  Define $\psi_n$ on $R_i$ by translating its lower left corner to that of $S_i$ and, relative to that common origin, separately rescaling the two coordinates.  Take $\psi$ to be the limit of this sequence of functions as $n$ increases.
Since each $\psi_n$ preserves probability by construction (each $R_i$ and $S_i$ comprises exactly $4^{-n}$ of the probability), so does $\psi.$  Thus, it converts one density into the other.  I leave the details of demonstrating convergence and the measurability of $\psi$ to you.  However, each $\psi_n$ is only piecewise linear: it potentially fails to be differentiable at the boundaries between the rectangles.

Reference
Roger B. Nelsen, An Introduction to Copulas, 2nd Ed. (Springer 2006), Sections 2.1 - 2.4.
