Computes the product moments of the truncated multivariate normal distribution, particularly for cases involving patterned variance-covariance matrices. It also has the capability to calculate these moments with arbitrary positive-definite matrices, although performance may degrade for high-dimensional variables.