This paper studies the impact of optimal sizing of photovoltaic distributed generators (PV-DGs) on a distribution system using different static load models (i.e., constant power, constant current, and constant impedance) and various power factor (PF) operations. A probabilistic approach with Monte Carlo simulation is proposed to obtain the optimal size of PV-DG. Monte Carlo simulation is applied to predict the solar radiation, ambient temperatures, and load demands. The objective is to minimize average system real power losses, with the power quality constraints not exceeding the limits, i.e. voltage and total harmonic voltage distortion (THDv) at the point of common coupling (PCC). A modified Newton method and a classical harmonic flow method are employed to calculate the power flow and THDv values, respectively. An actual 51-bus, medium-voltage distribution system in Thailand is employed as a test case. Results demonstrate that the proposed method performs well to provide the optimal size of PV-DG based on technical constraints. Further, the results show that the three static load models do not affect the optimal PV-DG size but the model has a different impact for various PF operations. PV-DGs may improve the voltage regulation and decrease the losses in distribution systems practically, but the THDv values could increase.