The large-scale integration of distributed energy resources significantly increases the uncertainty of harmonic power flow (HPF) in power systems. The traditional HPF calculation methods struggle to accurately capture the tail characteristics of harmonic distributions, thereby limiting the accuracy of HPF evaluation. Aiming at this problem, a probabilistic harmonic power flow (PHPF) calculation method for distribution networks based on improved point estimate method (PEM) is proposed. First, the traditional power flow method is used to calculate the system sensitivity discrimination index to identify the key harmonic sources. Second, using the hybrid sampling strategy and based on the Latin hypercube sampling, importance sampling is adopted specifically for the tail regions of critical harmonic sources. Finally, deterministic HPF calculations are performed on all samples, and the probability distributions of the outputs are characterized using kernel density estimation. Simulation results show that the mean error of the harmonic voltage amplitude of the proposed method is reduced by approximately 4.5% compared to the traditional PEM, the absolute value of the variance error is reduced by about 30%, while improving the fitting accuracy of the top 10% tail distribution. These results verify the effectiveness of the proposed method.