The expectation–maximization (EM) algorithm is a seminal method to calculate the maximum likelihood estimators (MLEs) for incomplete data. However, one drawback of this algorithm is that the asymptotic variance–covariance matrix of the MLE is not automatically produced. Although there are several methods proposed to resolve this drawback, limitations exist for these methods. In this paper, we propose an innovative interpolation procedure to directly estimate the asymptotic variance–covariance matrix of the MLE obtained by the EM algorithm. Specifically we make use of the cubic spline interpolation to approximate the first-order and the second-order derivative functions in the Jacobian and Hessian matrices from the EM algorithm. It does not require iterative procedures as in other previously proposed numerical methods, so it is computationally efficient and direct. We derive the truncation error bounds of the functions theoretically and show that the truncation error diminishes to zero as the mesh size approaches zero. The optimal mesh size is derived as well by minimizing the global error. The accuracy and the complexity of the novel method is compared with those of the well-known SEM method. Two numerical examples and a real data are used to illustrate the accuracy and stability of this novel method.