I am using Gram-Schmidt orthogonalization to create orthogonal polynomials with respect to an arbitrary weight function in a n-dimensional space (not always orthotope).
$P_i(\mathbf{x})=B_i(\mathbf{x})-\displaystyle\sum_{j=0}^{i-1}\frac{\left\langle B_i,P_j\right\rangle}{\left\langle P_j,P_j\right\rangle}P_j(\mathbf{x})$
Where $P_i$ is the i-th polynomial of my orthogonal basis and $\{B_i\}$ is an initial (and known) set of polynomials. My problem is that I have an inner product $\langle\cdot\rangle$ i.e. an integral in nD to perform numerically, the use of a quadrature rule requiters an orthogonal polynomial basis that I don't have yet since I am constructing it.
How do people generally perform this integral ? I was thinking about a Monte-Carlo integration but I may generate a lot of errors, especially in unbounded integration domains.