We present a method for reconstructing the photon number distribution from the homodyne statistics based on maximization of the likelihood function derived from the exact statistical description of a homodyne experiment. This method incorporates in a natural way the physical constraints on the reconstructed quantities, and the compensation for the nonunit detection efficiency.