How are the numbers of modes of marginal and joint distributions related? Can a unimodal multivariate distribution have a multimodal marginal distribution? 
If all marginal distributions are unimodal, can the multivariate distribution be multimodal?
 A: Yes and yes, and examples are easy to devise.
Let $(X,Y)$ have discrete distribution that has


*

*probability mass of $0.12$ at $(0,0)$ and $0.08$ at the $4$ points $(\pm 1,\pm 1)$.

*probability mass of $0.07$ at each of $(\pm 2, \pm 2)$ and $(\pm 3, \pm 3)$.
This is a unimodal bivariate distribution with a mode of $0.12$ at $(0,0)$
However, the marginal distributions 
are
$$\begin{array}{cccccccc}
-3& -2& -1& 0 & +1 & +2 & +3\\
0.14 & 0.14& 0.16 & 0.12 & 0.16 & 0.14 & 0.14
\end{array}$$
with twin peaks of $0.16$ at $\pm 1$ and least likely value $0$ where
the mode of the bivariate distribution is!

As an example of unimodal marginal distributions but a 
multimodal joint distribution, consider random variables $X$ and $Y$
that have identical unimodal marginal distributions 
$$\begin{array}{cccccccc}
-2& -1& 0 & +1 & +2 \\
0.05 & 0.2& 0.5 & 0.2 & 0.05
\end{array}$$
and bivariate distribution that has masses of $0.2$ at $(0,1), (1,0), (0,-1), (-1,0)$ and masses of $0.05$ at $(0,2), (2,0), (0,-2), (-2,0)$ which has 
$4$ modes.
A: One way to create a multimodal bivariate distribution out of unimodal marginals is with a mixture. The idea is that a mixture of Normal distributions can be unimodal, but bivariate versions of the same normals can have distinct peaks provided they are suitably correlated.
As a concrete example, let $f$ be an equal mixture of two bivariate Normal distributions.  Both have unit variances and the same correlation $\rho$, but let their means be $(a,-a)$ and $(-a,a)$ for some $a\ge 0$.  Provided $a\le 1,$ the marginals of $f$ (which are identical mixtures of two unit-variance Normals) will have just a single mode, but if we make $\rho$ close enough to $1$, $f$ will be distinctly bimodal. Here, for instance, are a contour and surface plot of $f$ for $a=1/2$ and $\rho=4/5$:

The common marginal has no hint of two modes:

The general technique illustrated here shows how a suitably chosen mixture of multivariate distributions can smooth out the marginals, causing them to be only unimodal. We see these kinds of things all the time when creating kernel density estimates of multivariate data: the density estimates have modes where the data cluster, but the marginals can be unimodal because the clusters do not align along any of the component directions.

To address the second question, observe that even a distribution all of whose marginals are Standard Normal can still be multimodal.
As a concrete example, let $f$ be the PDF of a bivariate Normal distribution with mean $(0,0)$, unit variances, and correlation $\rho$.  Define
$$g(x,y) = \begin{cases}
 f(-y,x) & -1<x<1\text{ and } -1<y<1 \\
 f(x,y)  & \text{otherwise}.
\end{cases}
$$
This rotates the graph of $g$ within the square $(-1,1)\times(-1,1)$ by 90 degrees. 
Because the marginals of $f$ are identical and symmetric around $0$, $g$ has the same marginals as $f$.  However, $g$ clearly has three modes--at $(0,0)$, $(1,1)$, and $(-1,-1)$--as this contour plot (with $\rho=4/5$) shows:

This example illustrates a general technique to answer questions about modes: cut pieces out of the PDF and move them around.
