SUMMARY: The fixed parameters of the nonlinear mixed effects model and the density of the random effects are estimated jointly by maximum likelihood. The density of the random effects is assumed to be smooth but is otherwise unrestricted. The method uses a series expansion that follows from the smoothness assumption to represent the density and quadrature to compute the likelihood. Standard algorithms are used for optimization. Empirical Bayes estimates of random coefficients are obtained by computing posterior modes. The method is applied to data from pharmacokinetics, and properties of the method are investigated by application to simulated data. © 1993 Biometrika Trust.