Bayesian nonlinear system identification for one of the major classes of dynamic model, the nonlinear autoregressive with exogenous input (NARX) model, has not been widely studied to date. Markov chain Monte Carlo (MCMC) methods have been developed, which tend to be accurate but can also be slow to converge. In this contribution, we present a novel, computationally efficient solution to sparse Bayesian identification of the NARX model using variational inference, which is orders of magnitude faster than the MCMC methods. A sparsity-inducing hyper-prior is used to solve the structure detection problem. Key results include: 1) successful demonstration of the method on low signal-to-noise ratio signals (down to 2 dB); 2) successful benchmarking in terms of speed and accuracy against a number of other algorithms: Bayesian LASSO, reversible jump MCMC, forward regression orthogonalization, LASSO, and simulation error minimization with pruning; 3) accurate identification of a real world system, an electroactive polymer; and 4) demonstration for the first time of numerically propagating the estimated nonlinear time-domain model parameter uncertainty into the frequency domain.