The statistical bias of numerically integrated statistical procedures

Maximum likelihood estimation of multivariate normal models and Bayesian posterior density functions can require calculation using numerical integration. Although the statistical theory is unchanged for nonanalytical integrals, the exorbitant cost of estimating such models frequently induces researchers to use a small number of integration points. If the number of points is insufficient to calculate essentially exact values of the function being maximized, the resulting estimators are asymptotically biased and inconsistent both for the values of the estimated parameters and for their variances. The forms of the biases are derived and described in this paper. The biases cannot be signed in general. The variance-covariance matrix is multiplied by an unestimable scalar, and the vector of true coefficients is added to an unestimable vector. When double integrals are involved, the number of points needed to calculate the integrals to desired accuracy is squared; with triple integrals, it is cubed, etc. The biases are additive, each integral contributing terms of similar form. An example of a multivariate normal model is presented. The paper shows that skimping on the high computational cost endangers the desirable statistical properties of MLE.