We present a new estimate of foreground emission in the Wilkinson Microwave Anisotropy Probe (WMAP) data, using a Markov chain Monte Carlo method. The new technique delivers maps of each foreground component for a variety of foreground models with estimates of the uncertainty of each foreground component, and it provides an overall goodness-of-fit estimate. The resulting foreground maps are in broad agreement with those from previous techniques used both within the collaboration and by other authors. We find that for WMAP data, a simple model with power-law synchrotron, free-free, and thermal dust components fits 90% of the sky with a reduced xi 2 Delta of 1.14. However, the model does not work well inside the Galactic plane. The addition of either synchrotron steepening or a modified spinning dust model improves the fit. This component may account for up to 14% of the total flux at the Ka band (33 GHz). We find no evidence for foreground contamination of the cosmic microwave background temperature map in the 85% of the sky used for cosmological analysis.