# Huber's M estimator for contaminated Gaussian noise

Huber discussed in this seminal paper "Robust Estimation of a Location Parameter" link that if we have some observations $x_i$ as follows: $$y_i = \theta + \nu_i, ~~i=1,\cdots,N, \tag{1}$$ where $\nu_i$ are given by the contaminated Gaussian distribution $$g(\nu) = (1-\epsilon) \phi(\nu) + \epsilon h(\nu). \tag{2}$$ Specifically, $\phi(\nu)$ is a Gaussian distribution, i.e., $\phi(\nu) = \frac{1}{\sqrt{2\pi}}e^{-\nu^2/2}$ and $h(\nu)$ can be any arbitrary symmetric distribution. Then, the M-estimator $\hat\theta$ of $\theta$ is given by the solution of $$\sum\limits_{i=1}^N \Psi(y_i - \theta) = 0, \tag{3}$$ where $$\Psi(x) = \left\{ \begin{array}{r@{\;\;}l} &x, ~~~|x|\leq k\\ &k , ~~~x>k \\ &-k, ~~~x<-k \end{array} \right. \tag{4}$$ and $k$ is a given number depending on $\epsilon$.

However, I could not find in his paper how exactly we can get $\hat\theta$ by solving (3). It is just mentioned in the paper that it can be solved by using some modified Winsorized mean. But I need the detailed steps. For example, if $k=1, y=[1,1.5,1,1.3]$, what is $\hat\theta$?

Any suggestion or reference will be appreciated.