In this paper, the probabilistic optimal power flow (POPF) problem is modeled as a black box uncertainty quantification problem. If distribution functions of POPF inputs are unknown, only samples are available, the multivariate Gaussian mixture model is suggested to recover the joint cumulative distribution function of POPF inputs, and a probability selection method is introduced to relate POPF inputs to an independent standard normal vector. Then, moment matching equations are proposed to characterize the uncertainty effects of POPF inputs on outputs, whereby some properties of the POPF problem are revealed, and a stratified quadrature rule (SQR) is proposed for POPF computation by combining a Hadamard matrix based quadrature rule and a discrete sine transformation matrix based quadrature rule. This newly developed SQR has a computational burden linear with the number of POPF inputs, and performs well in estimating means and standard deviations POPF outputs. Finally, case studies are conducted on a modified IEEE 118-bus system to illustrate the accuracy and efficiency of the proposed method.