This paper proposes a probabilistic optimal power flow (POPF) for unbalanced three-phase electrical distribution systems considering robust constraints. Wind velocity, solar irradiation, and demands are considered as exogenous random variables with dissimilar probability distribution functions. The proposed two-stage POPF model determines the optimal generation dispatch that minimizes the average production cost while satisfying probabilistically robust constraints. The robustness of the solution is adjusted through a parameter involving the first two moments of the random state variables. The $2m+1$ point estimate method (PEM) provides the scenarios for the proposed POPF model, which can be solved using commercial solvers. The model's accuracy is evaluated by comparing the results obtained with the PEM and with a classical scenario generator technique. The model is validated through two IEEE distribution feeders with 13 and 123 nodes, comprising dispatchable and renewable generation. Results show that the model can accurately estimate the mean and standard deviation of the state random variables when compared with results obtained through Monte Carlo simulations. Various levels of conservatism are assessed, demonstrating the model's flexibility and applicability.