For $M$ a simple surface, the non-linear statistical inverse problem of recovering a matrix field $Φ: M \to \mathfrak{so}(n)$ from discrete, noisy measurements of the $SO(n)$-valued scattering data $C_Φ$ of a solution of a matrix ODE is considered ($n\geq 2$). Injectivity of the map $Φ\mapsto C_Φ$ was established by [Paternain, Salo, Uhlmann; Geom.Funct.Anal. 2012]. A statistical algorithm for the solution of this inverse problem based on Gaussian process priors is proposed, and it is shown how it can be implemented by infinite-dimensional MCMC methods. It is further shown that as the number $N$ of measurements of point-evaluations of $C_Φ$ increases, the statistical error in the recovery of $Φ$ converges to zero in $L^2(M)$-distance at a rate that is algebraic in $1/N$, and approaches $1/\sqrt N$ for smooth matrix fields $Φ$. The proof relies, among other things, on a new stability estimate for the inverse map $C_Φ\to Φ$. Key applications of our results are discussed in the case $n=3$ to polarimetric neutron tomography, see [Desai et al., Nature Sc.Rep. 2018] and [Hilger et al., Nature Comm. 2018]