High-order discretizations and modern physics-informed and variational machine learning (ML) methods often require repeated numerical integration in multiple dimensions. Classical Gaussian-type quadrature rules are optimal in terms of exactness per point for polynomials, but their point locations are deterministic and highly structured, introducing systematic bias when the integrands depart from the target polynomial space (e.g., oscillatory or neural-network-based trial functions). This limits their robustness and increases overfitting when quadrature is embedded in training loops. We present a framework to construct stochastic quadrature rules with unbiasedness in expectation. A stochastic rule is defined as a set of deterministic quadratures equipped with probabilities. The method consists of two stages: First, we generate a large ensemble of randomized quadratures that remain exact on a polynomial space of degree p, with variability in point locations and weights. This variability is possible by increasing the number of points q in each rule, and thus expanding the degrees of freedom. This step is posed as an optimization problem and solved with neural-network-based parameterizations. Second, we unbias the ensemble by optimizing the probabilities so that the expected weighted distribution of points approximates the uniform (Lebesgue) measure on the integration domain. Practically, we minimize a differentiable histogram-based loss related to maximizing the entropy of the induced sampling distribution. We conducted numerical tests in 1D and 2D, with randomized rules exact up to degrees p=11 and p=5, respectively, on quadrilateral elements. While the initial randomized rules show biases approaching 20% for polynomials of degree 40, after probability optimization the bias is reduced to less than 1% in all cases. This method is applicable both to tensor-product and trunk-based spaces, can be adapted to other geometries, and has relevant practical implications for Deep-Ritz-type solvers.