Abstract Diffusion magnetic resonance imaging (dMRI) allows to estimate brain tissue microstructure as well as the connectivity of the white matter (known as tractography). Accurate estimation of the model parameters (by solving the inverse problem) is thus very important to infer the underlying biophysical tissue properties and fiber orientations. Although there has been extensive research on this topic with a myriad of dMRI models, most models use standard nonlinear optimization techniques and only provide an estimate of the model parameters without any information (quantification) about uncertainty in their estimation. Further, the effect of this uncertainty on the estimation of the derived dMRI microstructural measures downstream (e.g., fractional anisotropy) is often unknown and is rarely estimated. To address this issue, we first design a new deep-learning algorithm to identify the number of crossing fibers in each voxel. Then, at each voxel, we propose a robust likelihood-free deep learning method to estimate not only the mean estimate of the parameters of a multi-fiber dMRI model (e.g., the biexponential model), but also its full posterior distribution. The posterior distribution is then used to estimate the uncertainty in the model parameters as well as the derived measures. We perform several synthetic and in-vivo quantitative experiments to demonstrate the robustness of our approach for different noise levels and out-of-distribution test samples. Besides, our approach is computationally fast and requires an order of magnitude less time than standard nonlinear fitting techniques. The proposed method demonstrates much lower error (compared to existing methods) in estimating several metrics, including number of fibers in a voxel, fiber orientation, and tensor eigenvalues. The proposed methodology is quite general and can be used for the estimation of the parameters from any other dMRI model.