I think the simplest way to understand why we cannot have a joint probability distribution for two incompatible observables $A,B$ (meaning non-commuting) is the following. With a slight abuse of notation, the joint probability is defined as:
$$ P_{A,B}(a,b) := \mathrm{Prob}(A=a,B=b) $$
which means, it's the probability of $A$ having value $a$ and $B$ having value $b$. In quantum mechanics this means that there is an eigenstate of $A$ with eigenvalue $a$ which is at the same time an eigenstate of $B$ with eigenvalue $b$. But if $[A,B]\neq 0$ this is notoriously not possible.
Note that, conversely, if $A$ and $B$ commute it is always possible to find a common eigenbasis and so the above prescription works fine.
This argument, however does not answer the other part of the question. Which is:
Why can I not define a bona fide characteristic function (i.e. which is the Fourier transform of a probability density) in case of non-compatible observables? And perhaps, is there any way to amend this?
As pointed correctly by the OP this has to do with Bochner's theorem which tells precisely what requirements a characteristic function has to satisfy. I will state Bochner's theorem in the form needed for our purposes.
Theorem (Bochner) (Univariate case) $\chi (t)$ is the Fourier transform of a probability density $P(\omega)$ ($t,\omega \in\mathbb{R}$) if and only if for any $n$-tuple $t_1,t_2,\ldots t_n$ ($t_k \in \mathbb{R}, \ k=1,2,\ldots,n$) the $n\times n$ matrix with entries $\chi_{i,j}:= \chi(t_i-t_j)$ is non-negative definite (and hermitian).
Note: for the Multivariate $d$-dimensional generalization simply consider the obvious rephrasing with $t,\omega, t_k \in \mathbb{R}^d$.
A simple way to understand Bochner's theorem is the following. A matrix $\chi$ is non-negative definite if and only if it can be written as $\chi = A A^\dagger$.
Let $P(\omega)$ be the Fourier transform of $\chi$. Then
$$ \chi_{i,j} = \int d\omega e^{i(t_i-t_j) \omega} P(\omega) \ \ \ \ \ (0) $$
which we write as
\begin{align} \chi_{i,j} &= \int d\omega e^{it_i \omega} P(\omega) e^{-it_j \omega} \\ & = (A A^\dagger)_{i,j} \end{align}
with
$$ A_{i,\omega} := e^{it_i \omega} \sqrt{P(\omega)} $$
which we can do since $P(\omega)$ is non-negative. So $P(\omega)$ non-negative means that $\chi_{i,j}$ is a non-negative definite matrix. This characterizes characteristic functions in the classical case.
Let us now turn to the quantum mechanics and consider the multivariate case, i.e., we have several observables which I call $X_1, X_2, \ldots X_n$ with spectra in $\omega_1, \ldots, \omega_n$. The conjugate variables being $t_1,\ldots,t_n$, and the notation
$$tX:=\sum_{k=1}^n t_k X_k$$
In this case the wannabe characteristic function is
$$ \chi(t):=\mathsf{E}[ e^{itX} ]= \operatorname{Tr} ( e^{itX} \rho ) $$
for some quantum state $\rho$ (a normalized non-negative matrix). We want to check under which conditions $\chi_{i,j}=\chi(t_1-t_j)$ is non-negative definite as a matrix.
If the $X_k$ were mutually commuting operators we would have
$$ e^{i(t_i-t_j) X} = e^{it_i X} e^{-it_j X} \ \ \ \ \ \ (1) $$
and then we could write
\begin{align} \chi_{i,j} &=\operatorname{Tr} \left ( e^{i(t_i-t_j) X} \rho \right) \\ &=\operatorname{Tr} \left ( e^{it_i X} e^{-it_j X} \rho \right ) \\ &=\operatorname{Tr} \left (e^{-it_j X} \sqrt{\rho} \sqrt{\rho} e^{it_i X} \right ) \end{align}
Now define the matrix $A_{j,lq} := \left ( e^{-it_j X} \sqrt{\rho} \right )_{l,q} =: A_{i,\xi}$ where $\xi=(l,q)$. We have
$$ \chi_{i,j} = \sum_{lq} A_{j,lq} \overline{A_{j,lq}} $$
which is of the form $BB^\dagger$ and proves that $\chi_{i,j}$ is non-negative definite. Obviously this whole construction breaks down if $X_k$ are not mutually commuting.
Strictly speaking (as pointed out correctly by @AcuriousMind) this does not prove that the matrix $\chi_{i,j}$ is not non-negative definite for non-mutually commuting observables. For that one should find a counterexample, i.e. show that $\chi$ has a negative eigenvalue. However it does show where the argument breaks down.
Added edit
Here I present a counterexample for a single qubit. It can be shown that for $n=2$ the $\chi_{i,j}$ matrix is always non-negative definite. So to look for the first counterexample we must take $n=3$.
Consider the following problem with incompatible (non-commuting) observables given by $\sigma^x$ and $\sigma^z$. As for the state we pick $\rho = | 0\rangle\langle 0|$. Hence the putative characteristic function is
\begin{align} \chi(t_x,t_z) &:= \langle 0| e^{i (t_x \sigma^x +t_z \sigma^z)}| 0 \rangle \\ & = \cos \left (\sqrt{t_x^2 +t_z^2}\right ) + i\frac{t_z}{\sqrt{t_x^2 +t_z^2}} \sin \left ( \sqrt{t_x^2 +t_z^2}\right ) . \end{align}
Now form the matrix $\chi_{i,j}$ for $n=3$. It can be shown that $\chi$ has the form
$$ \chi = 1\!\mathrm{l} + \Gamma $$
where the matrix $\Gamma$ is hermitian and has zero on the diagonal. Since $\Gamma$ is traceless $\chi$ fails to be non-negative definite if $\Gamma$ has an eigenvalue smaller than $-1$.
For simplicity let's call $a_{ij} = t_x^i-t_x^j$ and $b_{ij} = t_z^i-t_z^j$. Now simply pick random $a_{ij}, b_{ij}$:
\begin{align} a_{12} & = 1 \ \ b_{12} = 0.5 \\ a_{13} & = 2 \ \ b_{13} = 1.3 \\ a_{23} & = 0.4 \ \ b_{23} = 0.9 \\ \end{align}
The eigenvalues of $\Gamma$ turn out to be $\{1.499, \ -1.221, \ -0.278\}$, which implies $\chi$ is not non-negative definite. This shows that the Fourier transform of $\chi(t_x,t_z)$ is not a (joint) probability distribution. $\square$