Robust estimators for count data I am looking for robust estimators for parameters of a processes producing count data:
$$
(n_1,...,n_K), n_i\in\mathbb{N}
$$
that is the underlying distribution is something like Poisson or Negative binomial, but possibly contaminated - with some samples exhibiting an excessively low or excessively high counts.
I have in mind something like an estimator along the lines of truncated mean or Winsorized mean (that is discarding some data), or an M-estimator along the lines of Huber loss giving lower weight to outliers.
Remarks

*

*Thread Robust estimation of Poisson distribution discusses solving a similar problem in R (robust package description). I am however more interested in the underlying math (which I intend to implement myself in Python.)

*Other R packages: poissonMT, robeth, robustbase

If the issue merely boils down to very high or very low observations, one would be tempted to just use a trimmed mean. The problem with that of course, is that your estimate may be biased. You could say bias shmias! We often trade an unbiased estimator for a biased one that achieves much better variance when the overall MSE is superior. In my toy example:
set.seed(123)
out <- replicate(1e3, {
  y <- rpois(1e3, 20)
  y[1:50] <- rpois(50, 5)
  y[51:100] <- rpois(50, 50)
  q <- quantile(y, c(0.025, 0.975))
  ysub <- y[y > q[1] & y < q[2]]
  c(mean(y), mean(ysub))
})

The untransformed mean outperforms the trimmed mean though the trimmed mean does have a lower bias
> rowMeans(out)
[1] 20.74856 20.32148
> rowMeans((out-20)^2)
[1] 0.5816081 0.1286018

If an unbiased estimate is truly what you're after, a neat trick could be to trim the data, and then estimate the resulting Poisson density as a truncated one using maximum likelihood. Consider the example below:
set.seed(123)
y <- rpois(1e4, 20)
 
## only 1% of the data formally outlying
y[1:50] <- rpois(50, 5)
y[51:100] <- rpois(50, 50)

## trim 5% and treat it like a truncated Poisson
q <- quantile(y, c(0.025, 0.975))
ysub <- y[y > q[1] & y < q[2]]


negloglik <- function(lambda, x, q) {
  -sum(dpois(x=x, lambda=lambda, log=T) - 
    log(ppois(q=q[2], lambda=lambda) - ppois(q=q[1], lambda=lambda))
  )
}

nlm(negloglik, 10, x=ysub, q=q)

