This paper presents a technique to obtain the probability density functions of network harmonic impedances in an unbalanced three-phase distribution system. In these systems the impedances can significantly vary as a function of time due to the presence of single and multi-phase capacitors, generators and filters which can be switched on and off and also due to the variations in the balanced and unbalanced loads. One of the objectives of the paper is to demonstrate the need to take into account the probabilistic nature of the system harmonic impedances. This facilitates the identification of the critical operational problems such as the resonance conditions which may go undetected when the same unbalanced systems are analysed under strictly deterministic scenarios. Numerical examples on the three phase unbalanced IEEE 13-bus distribution test systems are provided to illustrate the application of the proposed method.