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

Yes, there are.  To name just one, which I've had good experiences with, you can minimize the Cramer-von Mises distance between the empirical distribution and the theoretical distribution with estimated parameters, possibly throwing out a specified percentage of the lowest and highest data points (or not; see below.)  The C-vM distance is:
$$\omega^2 = \int_{-\infty}^{\infty} [F_n(x) - F(x;\theta)]^2\,\mathrm{d}F(x;\theta)$$
Obviously, in the case of a discrete distribution, you would use the sum, and you would integrate / sum over a prespecified subset of the data, making the appropriate adjustments in the calculation of $F(x;\theta)$.
Some general links are below:
The Consistency and Robustness of Modified Cramer-Von Mises and Kolmogorov-Cramer Estimators
Robust Weighted Cramer-von Mises estimators of location...
Robust estimation vis minimum distance methods
Robust Estimation with Exponentially Tilted Hellinger Distance
To see why this estimator is robust against small amounts of contamination (we will address the larger amounts of contamination case below), consider the function itself.  The range of values of the square term is $[0,1]$; no matter how big an outlier in a sample of size $N$ we see, the influence of it on the squared term is limited by the $1/N$ change it can inflict on the empirical CDF.  Likewise, as an outlier becomes more extreme, $\mathrm{d}F(x;\theta)$ goes to zero, further limiting its influence.  There still is some influence, as the outlier reduces the empirical CDF of all the non-outlier data points by a factor of $(N-1)/N$.
To illustrate this, we construct an example.  We generate 100 data points from a $\mathrm{Poisson(1)}$ distribution, and replace the $100^{th}$ with an increasing outlier.  We calculate the C-vM estimate of the parameter $\theta$, and plot how the estimate changes as the outlier increases:
cvm.calc <- function(theta, x) {
  ub <- max(qpois(0.9999, theta)+1, length(x))
  if (ub > length(x)) {
    x <- c(x, rep(1, ub-length(x)))
  }
  
  sum_over_range <- 0:(ub-1)
  sum((x-ppois(sum_over_range, theta))^2 * dpois(sum_over_range, theta))
}

x <- rpois(100, 1)
estimate <- rep(0,51)
for (outlier in 0:50) {
  x[100] <- outlier
  xtbl <- table(factor(x, 0:max(x)))
  x_edf <- cumsum(xtbl) / sum(xtbl)
  estimate[outlier+1] <- optimize(cvm.calc, interval=c(0.1,5), x=x_edf)$minimum
}

plot(estimate~I(0:50), type="b", pch=16,
     main="Sensitivity Curve - 1% Contaminated Poisson(1)",
     xlab="Outlier (1 % contamination)", 
     ylab="C-vM Estimate")

mean(x[1:99])  # The sample mean of the first 99 observations
[1] 1.030303
mean(x)        # The sample mean with the 100th observation = 50
[1] 1.52

with the result:

A similar plot, this time with 10% of the data replaced by outliers, results in:

with a sample mean in the outlier <- 50 case of 5.9.
In the more general case, where we expect larger amounts of contamination, we can always use a truncated version of both the empirical distribution function and the theoretical distribution function, e.g., dropping the top $\alpha \%$ of the data and adjusting the theoretical distribution accordingly.
An obvious question is: how efficient is it in the zero-outlier case?  A simulation will help to answer that in the $\mathrm{Poisson(1)}$ case:
cvm_estimate <- rep(0,10000)
for (i in seq_along(cvm_estimate)) {
  x <- rpois(100, 1)
  xtbl <- table(factor(x, 0:max(x)))
  x_edf <- cumsum(xtbl) / sum(xtbl)
  cvm_estimate[i] <- optimize(cvm.calc, interval=c(0.1,5), x=x_edf)$minimum
} 

mean((cvm_estimate - 1)^2)
[1] 0.01102528

which is only slightly worse than the MSE of the sample mean $(0.01)$.
