  The exploration of the low acceleration $a<a_{0}$ regime, where $a_{0}=1.2
\times 10^{-10}$m s$^{-2}$ is the acceleration scale of MOND around which
gravitational anomalies at galactic scale appear, has recently been extended to
the much smaller mass and length scales of local wide binaries thanks to the
availability of the {\it Gaia} catalogue. Statistical methods to test the
underlying structure of gravity using large samples of such binary stars and
dealing with the necessary presence of kinematic contaminants in such samples
have also been presented. However, an alternative approach using binary samples
carefully selected to avoid any such contaminants, and consequently much
smaller samples, has been lacking a formal statistical development. In the
interest of having independent high quality checks on the results of wide
binary gravity tests, we here develop a formal statistical framework for
treating small, clean, wide binary samples in the context of testing
modifications to gravity of the form $G \to \gamma G$. The method is validated
through extensive tests with synthetic data samples, and applied to recent {\it
Gaia} DR3 binary star observational samples of relative velocities and internal
separations on the plane of the sky, $v_{2D}$ and $r_{2D}$, respectively. Our
final results for a high acceleration $r_{2D}<0.01$pc region are of
$\gamma=1.000 \pm 0.096$, in full accordance with Newtonian expectations. For a
low acceleration $r_{2D}>0.01$pc region however, we obtain $\gamma=1.5 \pm
0.2$, inconsistent with the Newtonian value of $\gamma=1$ at a $2.6 \sigma$
level, and much more indicative of MOND AQUAL predictions of close to
$\gamma=1.4$.
