  We use a separable mode expansion estimator with WMAP data to estimate the
bispectrum for all the primary families of non-Gaussian models. We review the
late-time mode expansion estimator methodology which can be applied to any
non-separable primordial and CMB bispectrum model, and we demonstrate how the
method can be used to reconstruct the CMB bispectrum from an observational map.
We extend the previous validation of the general estimator using local map
simulations. We apply the estimator to the coadded WMAP 5-year data,
reconstructing the WMAP bispectrum using $l<500$ multipoles and $n=31$
orthonormal 3D eigenmodes. We constrain all popular nearly scale-invariant
models, ensuring that the theoretical bispectrum is well-described by a
convergent mode expansion. Constraints from the local model $ \fnl=54.4\pm
29.4$ and the equilateral model $\fnl=143.5\pm 151.2$ ($\Fnl = 25.1\pm 26.4$)
are consistent with previously published results. (Here, we use a nonlinearity
parameter $\Fnl$ normalised to the local case, to allow more direct comparison
between different models.) Notable new constraints from our method include
those for the constant model $\Fnl = 35.1 \pm 27.4 $, the flattened model $\Fnl
= 35.4\pm 29.2$, and warm inflation $\Fnl = 10.3\pm 27.2$. We investigate
feature models surveying a wide parameter range in both the scale and phase,
and we find no significant evidence of non-Gaussianity in the models surveyed.
We propose a measure $\barFnl$ for the total integrated bispectrum and find
that the measured value is consistent with the null hypothesis that CMB
anisotropies obey Gaussian statistics. We argue that this general bispectrum
survey with the WMAP data represents the best evidence for Gaussianity to date
and we discuss future prospects, notably from the Planck satellite.
