Due to the ever-increasing use of non-linear loads and their undesired effects on distribution systems operation, harmonic analysis should be taken into consideration. On the other hand, the probabilistic nature of power systems makes it necessary to consider the harmonic analysis in a probabilistic environment. In this paper, a data clustering based algorithm is used for probabilistic assessment of harmonic load flow, for the first time. Despite the previous probabilistic harmonic load flow (PHLF), in which uncertainties are considered on grid connected renewable generations, load demands, generators, transmission lines probable failure, etc., this paper considers uncertainties on the location and non-linear load portion of nodal loads. Moreover, an organized PHLF algorithm is formulated in this paper. In order to show the superior abilities of the proposed method, the method is applied on the IEEE 37 node test system and the results are compared by the Monte Carlo simulation (MCS) method.