  We introduce qustochastic matrices as the bistochastic matrices arising from
quaternionic unitary matrices by replacing each entry with the square of its
norm. This is the quaternionic analogue of the unistochastic matrices studied
by physicists. We also introduce quaternionic Hadamard matrices and
quaternionic mutually unbiased bases (MUB). In particular we show that the
number of MUB in an n-dimensional quaternionic Hilbert space is at most 2n+1.
The bound is attained for n=2. We also determine all quaternionic Hadamard
matrices of size at most 4.
