In pre-stack seismic amplitude variation with incident angle (AVA) inversion, observational noise is typically modeled as a Gaussian distribution with 0 mean and known (constant) covariance. Commonly, the forward equation of AVA inversion is based on Zoeppritz equations. But in practice, due to the complexity of the exact Zoeppritz equations, approximations are frequently used in industry to rapidly calculate AVA inversion. These approximations assume small elastic parameter contrasts. There is no significant loss of precision when the assumption is approximately satisfied. The approximations become less accurate for large elastic parameter contrasts. The essential reason is that the approximations can introduce model errors. This study employs a Bayesian approach to comprehensively characterize the statistics of model errors between the exact Zoeppritz equations and their approximations. Specifically, the Bayesian approximation error method is used to handle model errors, compensating for errors introduced by approximations of the Zoeppritz equations. Then, these model-error statistics are applied to AVA inversion for elastic parameters based on Bayesian theory. The synthetic data tests and real seismic data applications showed that the proposed AVA inversion method can improve the precision and accuracy of inverted elastic parameters to a certain extent, with little additional computational cost compared to AVA inversion based on approximations without model-error compensation.