  Standard confidence intervals employed in applied statistical analysis are
usually based on asymptotic approximations. Such approximations can be
considerably inaccurate in small and moderate sized samples. We derive accurate
confidence intervals based on higher-order approximate quantiles of the score
function. The coverage approximation error is $O(n^{-3/2})$ while the
approximation error of confidence intervals based on the asymptotic normality
of MLEs is $O(n^{-1/2})$. Monte Carlo simulations confirm the theoretical
findings. An implementation for regression models and real data applications
are provided.
