Motivation for Wigner phase space distribution Most sources say that Wigner distribution acts like a joint phase-space distribution in quantum mechanics and this is justified by the formula 
$$\int_{\mathbb{R}^6}w(x,p)a(x,p)dxdp= \langle \psi|\hat{a}\psi\rangle$$
as this allows us to compute the expectation value for the operator $\hat{a}$ corresponding to the physical quantity $a$. However, without the knowledge of this formula, how did Wigner come up with this definition:
$$w(x,p)=\frac{1}{(2\pi)^3}\int_{\mathbb{R}^3}\psi\left(x-\frac{v}{2}\right)\psi^*\left(x+\frac{v}{2}\right)e^{iv \cdot p}dv \quad ?$$
I would be greatly indebted for any mathematical motivation that anyone could provide. Similarly, I also wonder at the mathematical motivation for the Weyl quantization.
 A: I) Let us for simplicity work in 1D with $\hbar=1$. (The generalization to higher dimensions is straightforward.) Moreover, let us for simplicity take an operator $\hat{f}(\hat{X},\hat{P})$ without any ordering ambiguities, i.e., each monomial term in the symbol $f(x,p)$ depends only on either $x$ or $p$, but not on both. Then one possible motivation of Wigner's phase space distribution 
$$
\tag{1} w(x,p)~:=~ \int_{\mathbb{R}}\!  {dy\over2\pi}e^{ipy}\psi^{*}(x+\frac{y}{2}) \psi(x-\frac{y}{2})
$$
goes as follows. The expectation value of the operator $\hat{f}(\hat{X},\hat{P})$ in the Schrödinger position representation 
$$\tag{2}  \hat{X} ~\longrightarrow~x, \qquad  \hat{P} ~\longrightarrow~  -i \frac{\partial}{\partial x}, $$
reads 
$$\langle\psi| \hat{f}(\hat{X},\hat{P})|\psi \rangle ~\stackrel{(2)}{=}~ \int_{\mathbb{R}} \!dx~ \psi^{*}(x) f\left(x, -i \frac{\partial}{\partial x}\right) \psi(x) \qquad\qquad  $$
$$ ~\stackrel{\begin{matrix}\text{substitute}\\ x=x^{\pm}\end{matrix}}{=}~ \int_{\mathbb{R}^{2}}\! dx^{+}dx^{-}~ \delta(x^{+}-x^{-})\psi^{*}(x^{+}) f\left(\frac{x^{+}+x^{-}}{2}, -i \frac{\partial}{\partial x^{-}}\right) \psi(x^{-}) $$
$$ ~\stackrel{\delta\text{-fct}}{=}~ \int_{\mathbb{R}^{3}}\!  {dx^{+}dx^{-}dp\over2\pi}~ e^{ip(x^{+}-x^{-})}\psi^{*}(x^{+}) f\left(\frac{x^{+}+x^{-}}{2}, -i \frac{\partial}{\partial x^{-}}\right) \psi(x^{-}) $$
$$ ~\stackrel{\text{int. by parts}}{=}~\int_{\mathbb{R}^{3}}\!  {dx^{+}dx^{-}dp\over2\pi}~ e^{ip(x^{+}-x^{-})}\psi^{*}(x^{+}) f\left(\frac{x^{+}+x^{-}}{2}, p\right) \psi(x^{-}) $$
$$ ~\stackrel{\begin{matrix}\text{substitute}\\ x^{\pm} = x\pm \frac{y}{2}\end{matrix}}{=}~\int_{\mathbb{R}^{3}}\!  {dx~dy~dp\over2\pi}~ e^{ipy}\psi^{*}(x+\frac{y}{2}) f(x, p) \psi(x-\frac{y}{2}) $$
$$ \tag{3}  ~\stackrel{(1)}{=}~\int_{\mathbb{R}^{2}}\!  dx~dp~w(x,p) f(x, p)  .$$
That's the motivation! 
II) For more general operators $\hat{f}(\hat{X},\hat{P})$,
we leave it for OP to show that if $f(x,p)$ is interpreted as the Weyl-symbol of the operator $\hat{f}(\hat{X},\hat{P})$, then the equation
$$\tag{3'} \langle\psi| \hat{f}(\hat{X},\hat{P}) |\psi\rangle ~=~\int_{\mathbb{R}^{2}}\!  dx~dp~w(x,p) f(x, p)$$
continues to hold [at least for a sufficiently well-behaved function $f(x, p)$].
III) One important virtue, from a physics perspective, of the Weyl-ordering (as opposed to other ordering prescriptions) is that the operator $\hat{f}(\hat{X},\hat{P})$ formally becomes Hermitian for real functions $f:\mathbb{R}^2 \to\mathbb{R}$ and two Hermitian operators $\hat{X}$ and $\hat{P}$. Recall that Hermitian operators correspond to physical observables in quantum mechanics. For Weyl-ordering, see also e.g. this Phys.SE post. 
A: I want to add an interesting detail to the excellent answer by Qmechanic which I recently learned from Cosmas Zachos, here.
Royer showed in his 1977 paper, that the Wigner function $w$ is given by the expectation value of the parity operator $\Pi_{\mathbf{r}\mathbf{p}}$ which inverts the phase space around the point $\mathbf{z}=(\mathbf{r},\mathbf{p})$, i.e., he showed that
$$
w(\mathbf{r},\mathbf{p}) \propto \langle \psi| \Pi_{\mathbf{r}\mathbf{p}}| \psi \rangle ,
$$
with (in one dimension for simplicity)
$$
\begin{align}
\Pi_{rp} &= \int\text{d} s~ e^{-2ips/\hbar} ~| r-s\rangle \langle r+s |
\\ 
&=
\int\text{d} k~ e^{-2ikr/\hbar} ~| p+k\rangle \langle p-k |,
\end{align}
$$
where $|r\rangle$ and $|p\rangle$ are eigenstates of position $\hat{r}$ and momentum $\hat{p}$. He proved that this operator has the property
$$
\begin{align}
\Pi_{rp} (\hat{r} -r) \Pi_{rp} & = -(\hat{r} -r) \\
\Pi_{rp} (\hat{p} -p) \Pi_{rp} & = -(\hat{p} -p) .
\end{align}
$$
As Royer states himself:

[The Wigner function] is proportional to the overlap of $\psi$ with its mirror image about $r,p$, which is clearly a measure of how much $\psi$ is "centered" around $r,p$.

Although this is not an answer to the question as Wigner himself probably did not know about this, it still adds some physical depth to the problem and could serve as a modern motivation. 
Something which you can immediately understand based on this interpretation is how the Wigner function is bounded and what these bounds represent. Royer showed that if $\psi$ is antisymmetric  about $r,p$, the Wigner function attains its lower bound.
