We propose a method to obtain the maximum likelihood (ML) parameter estimation of the Gamma-Gamma (Γ-Γ) distribution representing the free space optical (FSO) channel irradiance fluctuations. The proposed method is based on the expectation maximization (EM) algorithm and the generalized Newton method using a non-quadratic approximation. The numerical results show that, for all turbulence conditions, the proposed ML method is more accurate than the fractional moments (FMOM) method and the numerical ML method (two dimensional numerical maximization of log-likelihood function using the Nelder-Mead algorithm). Moreover, the proposed ML is a fast and stable iterative method, because the iterations always converge to the global optimum with high convergence rate.