We introduce a new method (to be referred to as the variational method, VM) for estimating the parameters of Gibbs distributions with random variables ("spins") taking values in a Euclidean space $\mathbb{R}^n, n \geq 1$, from complete or degraded data. The method applies also to the case of iid random variables governed by exponential families, and appears to be new even in this case. For complete data, the VM is computationally more efficient than, and as reliable as, the maximum pseudo-likelihood method. For incomplete data, the VM leads to an estimation procedure reminiscent of, but simpler than, the EM algorithm. In the former case, we show that under natural assumptions a certain form of the variational estimators is strongly consistent and asymptotically normal. We also present numerical experiments that demonstrate the computational efficiency and accuracy of the variational estimators.