This paper presents a method to forecast the probability distribution function (PDF) of the generated power of PV systems based on the higher order Markov chain (HMC). Since the output power of the PV system is highly influenced by ambient temperature and solar irradiance, they are used as important features to classify different operating conditions of the PV system. The classification procedure is carried out by applying the pattern discovery method on the historical data of the mentioned variables. An HMC is developed based on the categorized historical data of PV power in each operating point. The 15-min ahead PDF of the PV output power is forecasted through the Gaussian mixture method (GMM) by combining several distribution functions and by using the coefficients defined based on parameters of the HMC-based model. In order to verify the proposed method, the genetic algorithm is applied to minimize a well-defined objective function to achieve the optimal GMM coefficients. Numerical tests using real data demonstrate that the forecast results follow the real probability distribution of the PV power well under different weather conditions.