  Symplectic integrators are widely used for the study of planetary dynamics
and other $N$-body problems. In a study of the outer Solar system, we
demonstrate that individual symplectic integrations can yield biased errors in
the semi-major axes and possibly other orbital elements. The bias is resolved
by studying an ensemble of initial conditions of the outer Solar system. Such
statistical sampling could significantly improve measurement of planetary
system properties like their secular frequencies. We also compared the
distributions of action-like variables between high and low accuracy
integrations; traditional statistical metrics are unable to distinguish the
distribution functions.
