Calculating bivariate orthonormal polynomials by recurrence

Emerson gave recurrence formulae for the calculation of orthonormal polynomials for univariate discrete random variables. He claimed that as these were based on the Christoffel–Darboux recurrence relation they were more efficient than those based on the Gram–Schmidt method. This approach was generalised by Rayner and colleagues to arbitrary univariate random variables. The only constraint was that the expectations needed are well-defined. Here the approach is extended to arbitrary bivariate random variables for which the expectations needed are well-defined. The extension to multivariate random variables is clear.