A method for moving least squares interpolation and differentiation is presented in the framework of orthogonal polynomials on discrete points. This yields a robust and efficient method which can avoid singularities and breakdowns in the moving least squares method caused by particular configurations of nodes in the system. The method is tested by applying it to the estimation of first and second derivatives of test functions on random point distributions in two and three dimensions and by examining in detail the evaluation of second derivatives on one selected configuration. The accuracy and convergence of the method are examined with respect to length scale (point separation) and the number of points used. The method is found to be robust, accurate, and convergent.