number (2200) of mock catalogues allows a robust determination of the error covariance function between the bulk flow estimates at different radii. Since  is significantly superior to MLE, we restrict the comparison with models to  reconstruction. We use all components of $\vace$ estimated at 8 values of $r$ ranging from $r=30\hmpc$ to $100\hmpc$ in steps of $10\hmpc$. The reason for not considering smaller radii is that the bulk is most robustly constrained independent of the assumed basis functions at $r>30\hmpc$. We denote the set of  reconstructed cartesian components at these 8 values of $r$ by $\cBt$ and the corresponding underlying true quantities by ${\cal B}_t$. We write the probability for observing the set $\cBo$ as $$\label{eq:pr}
 P( \cBo)=    \int \dd \cBt P(\cBo|\cBt) P(\cBt) \;$$ where the probability $P(\cBt) $ for the underlying $\cBt$ is computed within the framework of a cosmological model. Here, we adopt the $\Lambda$CDM model. For gaussian velocity fields, the calculation of $P(\cBt)$ is easily done by integrating standard analytic expressions involving the power spectrum. We assume that the probability $P(\cBo|\cBt)$ for $\cBo$ given $\cBt$ is gaussian with error covariance matrix computed from the 2200 mocks. Under these assumptions, the expression (\[eq:pr\]) yields $$\label{eq:pbo}
 P( \cBo)=\frac{1}{\sqrt{(2\pi)^d |\Sigma|}}{\rm 