Summary: We develop statistical methodology for a popular brain imaging technique HARDI based on the high order tensor model by E. Özarslan and T. Mareci [“Generalized diffusion tensor imaging and analytical relationships between diffusion tensor imaging and high angular resolution diffusion imaging”, Magn. Reson. Med. 50, No. 5, 955–965 (2003; doi:10.1002/mrm.10596)]. We investigate how uncertainty in the imaging procedure propagates through all levels of the model: signals, tensor fields, vector fields, and fibers. We construct asymptotically normal estimators of the integral curves or fibers which allow us to trace the fibers together with confidence ellipsoids. The procedure is computationally intense as it blends linear algebra concepts from high order tensors with asymptotical statistical analysis. The theoretical results are illustrated on simulated and real datasets.{ }This work generalizes the statistical methodology proposed for low angular resolution diffusion tensor imaging by O. Carmichael and L. Sakhanenko [“Estimation of integral curves from noisy diffusion tensor data”, Preprint], to several fibers per voxel. It is also a pioneering statistical work on tractography from HARDI data. It avoids all the typical limitations of the deterministic tractography methods and it delivers the same information as probabilistic tractography methods. Our method is computationally cheap and it provides well-founded mathematical and statistical framework where diverse functionals on fibers, directions and tensors can be studied in a systematic and rigorous way.