  As modern-day precision cosmology aims for statistical uncertainties of the
percent level or lower, it becomes increasingly important to reconsider
estimator assumptions at each step of the process, and their consequences on
the statistical variability of the scientific results.
  We compare $L^1$ regression statistics to the weighted mean, the canonical
$L^2$ method based on Gaussian assumptions, for inference of the weak
gravitational shear signal from a catalog of background ellipticity
measurements around a sample of clusters, in many recent analyses a standard
step in the process.
  We use the shape measurements of background sources around 6925 AMICO
clusters detected in the KiDS 3rd data release. We investigate the robustness
of our results and the dependence of uncertainties on the signal-to-noise
ratios of the background source detections. Using a halo model approach, we
derive lensing masses from the estimated excess surface density profiles.
  The highly significant shear signal allows us to study the scaling relation
between the $r$-band cluster luminosity $L_{200}$, and the derived lensing mass
$M_{200}$. We show the results of the scaling relations derived in 13 bins in
$L_{200}$, with a tightly constrained power law slope of $\sim 1.24\pm 0.08$.
We observe a small, but significant relative bias of a few percent in the
recovered excess surface density profiles between the two regression methods,
which translates to a $1\sigma$ difference in $M_{200}$. The efficiency of
$L^1$ is at least that of the weighted mean, relatively increasing with higher
signal-to-noise shape measurements.
  Our results indicate the relevance of optimizing the estimator for infering
the gravitational shear from a distribution of background ellipticities. The
interpretation of measured relative biases can be gauged by deeper
observations, while increased computation times remain feasible.
