Harmonic distortion in electrical power systems is responsible for several technical problems. Harmonic load flow (HLF) methods have been employed in order to predict and solve many of the problems from a deterministic point of view. Moreover, methods based on probability theory have been developed to deal with the uncertainties in power systems and the random nature of harmonics. Unfortunately, since many of these methods are based on linearized models or some simplifying assumptions, they do not have acceptable accuracy. This paper proposes the use of an efficient sampling method, median Latin hypercube sampling, combined with an improved kernel density estimator, in Monte Carlo simulation for probabilistic HLF calculation. The proposed method has been applied to the well-known IEEE 14-bus harmonic test system to evaluate the harmonic probability density functions of output random variables. The simulation results clearly show that it guarantees a reasonable execution time as well as acceptable accuracy.