High-quality nanomechanical oscillators can sensitively probe force, mass, or displacement in experiments bridging the gap between the classical and quantum domain. Dynamics of these stochastic systems is inherently determined by the interplay between acting external forces, viscous dissipation, and random driving by the thermal environment. The relevance of inertia then dictates that both position and momentum must, in principle, be known to fully describe the system, which makes its quantitative experimental characterization rather challenging. We introduce a general method of Bayesian inference of the force field and environmental parameters in stochastic inertial systems that operates solely on the time series of recorded noisy positions of the system. The method is first validated on simulated trajectories of model stochastic harmonic and anharmonic oscillators with damping. Subsequently, the method is applied to experimental trajectories of particles levitating in tailored optical fields and used to characterize the dynamics of particle motion in a nonlinear Duffing potential, a static or time-dependent double-well potential, and a nonconservative force field. The presented inference procedure does not make any simplifying assumptions about the nature or symmetry of the acting force field and provides robust results with trajectories 2 orders of magnitude shorter than those typically required by alternative inference schemes. In addition to being a powerful tool for quantitative data analysis, it can also guide experimentalists in choosing appropriate sampling frequency (at least 20 measured points per single characteristic period) and length of the measured trajectories (at least 10 periods) to estimate the force field and environmental characteristics with a desired accuracy and precision.