Precise radial velocity measurements have led to the discovery of ~100\nextrasolar planetary systems. We investigate the uncertainty in the orbital\nsolutions that have been fit to these observations. Understanding these\nuncertainties will become more and more important as the discovery space for\nextrasolar planets shifts to longer and longer periods. While detections of\nshort period planets can be rapidly refined, planets with long orbital periods\nwill require observations spanning decades to constrain the orbital parameters\nprecisely. Already in some cases, multiple distinct orbital solutions provide\nsimilarly good fits, particularly in multiple planet systems. We present a\nmethod for quantifying the uncertainties in orbital fits and addressing\nspecific questions directly from the observational data rather than relying on\nbest fit orbital solutions. This Markov chain Monte Carlo (MCMC) technique has\nthe advantage that it is well suited for the high dimensional parameter spaces\nnecessary for the multiple planet systems. We apply the MCMC technique to\nseveral extrasolar planetary systems, assessing the uncertainties in orbital\nelements for several systems. Our MCMC simulations demonstrate that for some\nsystems there are strong correlations between orbital parameters and/or\nsignificant non-Gaussianities in parameter distributions, even though the\nmeasurement errors are nearly Gaussian. Once these effects are considered the\nactual uncertainties in orbital elements can be significantly larger or smaller\nthan the published uncertainties. We also present simple applications of our\nmethods such as predicting the times of possible transits for GJ 876.\n