Single photon count distributions at the outputs of a Mach-Zehnder interferometer are highly dependent on the phase difference between the two arms of the apparatus, so photocount records can be used to estimate an unknown optical phase. Experimental results have shown that phase information can be inferred from the data even at low mean input intensities of ⪅10 photons. Here we consider the optimal Bayesian parameter estimate for the model used by a previous frequentist analysis. In the limit of high input intensity the observation is well-modelled as an intensity measurement of a classical wave in Gaussian measurement noise, however in the low intensity regime this simple model does not apply. It seems that no generally valid closed-form estimator for this problem exists without making the Gaussian approximation. Moreover the 2π
radian ambiguity in phase estimates presents a problem in defining cost functions to compare estimator performance, and we consider two reasonable alternatives. We present a numerical study of the experiment along with a novel numerical approach to this problem which out-performs the existing estimators and with a squared error that approaches the Cramer-Rao lower bound previously derived in the literature, in the limit where that bound is valid.