[en] The equilibrium quasiprobability density function W(θ, ψ) of spin orientations in a representation (phase) space of the polar and azimuthal angles (θ, ψ) (analogous to the Wigner distribution for translational motion of a particle) is given by a finite series of spherical harmonics in the spin number and their associated statistical moments so allowing one to calculate W(θ, ψ) for an arbitrary spin system in the equilibrium state described by the canonical distribution ρ-hat=e^-βH-hat_S/Tr(e^-βH-hat_S). The system with Hamiltonian H-hat_S=-γℎH.S-hat-BS-hat_Z² is treated as a particular example (γ is the gyromagnetic ratio, ℎ is Planck's constant, H represents an external magnetic field and B represents an internal field parameter). For a uniaxial system with H-hat_S=-γℎHS-hat_Z-BS-hat_Z², the solution may be given in the closed form