  We introduce and carefully study a natural probability measure over the
numerical range of a complex matrix $A \in M_n(\C)$. This numerical measure
$\mu_A$ can be defined as the law of the random variable $<AX,X> \in \C$ when
the vector $X \in \C^n$ is uniformly distributed on the unit sphere. If the
matrix $A$ is normal, we show that $\mu_A$ has a piecewise polynomial density
$f_A$, which can be identified with a multivariate $B$-spline. In the general
(nonnormal) case, we relate the Radon transform of $\mu_A$ to the spectrum of a
family of Hermitian matrices, and we deduce an explicit representation formula
for the numerical density which is appropriate for theoretical and
computational purposes. As an application, we show that the density $f_A$ is
polynomial in some regions of the complex plane which can be characterized
geometrically, and we recover some known results about lacunas of symmetric
hyperbolic systems in $2+1$ dimensions. Finally, we prove under general
assumptions that the numerical measure of a matrix $A \in M_n(\C)$ concentrates
to a Dirac mass as the size $n$ goes to infinity.
