This paper proposes a new probabilistic load flow (PLF) method that considers the non-Gaussianity and the correlation of nodal input variables. The load flow calculation plays an important role to evaluate power system conditions. In recent years, power systems become more complicated due to renewable energy such as wind power generators and PV systems. As a result, the impotence of PLF has been reevaluated as one of the promising techniques for handling the network uncertainties. This paper focuses on the applications of Monte Carlo Simulation (MCS) to PLF so that the effect of the correlations of input variables and the network uncertainties is evaluated through the nonlinear equation. In this paper, new techniques are introduced into MCS-PLF to evaluate more accurate solutions. The maximum likelihood estimation for a probability density function (PDF) of input variables is carried out to construct the non-Gaussian distribution model of input variables with the correlation by the Deterministic Annealing Expectation Maximization (DAEM) algorithm. Also, the Metropolis-Hastings sampling is used to generate random numbers that correspond to complicated multivariate probability density functions. The proposed method is successively applied to a system with wind power generators.