Probabilistic load flow (PLF) allows to evaluate uncertainties introduced by renewable energy sources on system operation. Ideally, the PLF calculation is implemented for an entire grid requiring all the parameters of the transmission lines and node load/generation to be available. However, in a multi-regional interconnected grid, the independent system operators (ISOs) across regions may not share the parameters of their respective areas with other ISOs. Consequently, the challenge is how to identify the functional relationship between the flows in the regional grid and the uncertain power injections of renewable generation sources across regions without full information about the entire grid. To overcome this challenge, we first propose a privacy-preserving distributed accelerated projection-based consensus algorithm for each ISO to calculate the corresponding coefficient matrix of the desired functional relationship. Then, we leverage a privacy-preserving accelerated average consensus algorithm to allow each ISO to obtain the corresponding constant vector of the same relationship. Using the two algorithms, we finally derive a privacy-preserving distributed PLF method for each ISO to analytically obtain its regional joint PLF in a fully distributed manner without revealing its parameters to other ISOs. The correctness, effectiveness, and efficiency of the proposed method are verified through a case study on the IEEE 118-bus system.