  We propose a method to define quasiprobability distributions for general
spin-$j$ systems of dimension $n=2j+1$, where $n$ is a prime or power of prime.
The method is based on a complete set of orthonormal commuting operators
related to Mutually Unbiased Bases which enable (i) a parameterisation of the
density matrix and (ii) construction of measurement operators that can be
physically realised. As a result we geometrically characterise the set of
states for which the quasiprobability distribution is non-negative, and can be
viewed as a joint distribution of classical random variables assuming values in
a finite set of outcomes. The set is an $(n^2-1)$-dimensional convex polytope
with $n+1$ vertices as the only pure states, $n^{n+1}$ number of higher
dimensional faces, and $n^3(n+1)/2$ edges.
