Load flow studies are crucial in investigations of operation and planning problems in the power systems. Traditional methods for determining load flow are based on deterministic approaches. However, the parameters of a power system (such as load and renewable power generation) may be uncertain. An exact probabilistic load flow (PLF) study requires a long CPU time due to many convolution computations involving probability density functions. This paper proposes a novel PLF method that is based on Taguchi's orthogonal arrays. The proposed method utilises a few deterministic load flow solutions that are obtained using Taguchi's method to estimate the means and standard deviations of bus voltages, phase angles, line flows, and other metrics. A load flow calculation corresponds to an experiment in Taguchi's method. An optimal experiment is also specified by considering the largest deviation from the nominal load flow solution. A 25-bus standalone power system and a modified Institute of Electrical and Electronics Engineers (IEEE) 118-bus system are tested. The simulation results show that the proposed method not only requires fewer deterministic load flow solutions to perform PLF analysis than the traditional point-estimate method but also yields accurate means and standard deviations of bus voltages and line flows.