My problem is a parameters' estimation of highly contaminated normal distribution.

Typical estimators of variance such as IQR or MAD are based on the assumption that there are not that much outliers. But what if we have a lot of outliers, probably, more than 50% of points, but all of them are quite distant from the center of normal distribution so they affect only the tails? Can we, looking at the central part (that is not affected by outliers), estimate parameters of the normal distribution?

It can be the case when the measuring device 1) outputs normally distributed values if it is functioning properly, 2) outputs quite big values weirdly distributed when broken (and it can be broken >50% of cases). I have such situation in biology so it is not just my interest, it has a practical meaning.

Additional assumptions: In my situation the peak of normal distribution of interest is well-defined and can always be found with the mode estimation algorithm, so there is no problem with estimation of the center. I can also find the points when the outliers start to appear (vertical lines on the attached plot). On the plot I showed the uniformly distributed outliers, however, there is no assumption on their distribution, so it can not be solved just with EM mixture algorithm. I should really estimate the parameters based on something like kernel density estimation and I know that I have to cut tails since they are spoiled.

enter image description here UPDATE: reformulated question:

enter image description here

I draw 6 points on the plot of standard normal distribution (ABCDEF). The distance between A and C is an IQR of the distribution. We can infer standard deviation from the length of $|AC|$ - we should just divide $|AC|$ with the ~1.3489. The question now is: how can we infer a standard deviation, if we move the perpendiculars $AB$ and $DC$ on the $x$-axis? We will know coordinates of all 6 points on both $y$ and $x$ axis, and will know that $EF$ is in the center of the normal distribution.

Draft of the solution: $x$ is the value we are looking for (standard deviation). $dnorm$ is function that returns likelihood. $\hat \mu$ can be estimated from the data. When I use letter $A$, I mean the $x$-coord of the point.

$\frac{|AB|}{|EF|} = \frac{dnorm(A, \mu=\hat \mu, \sigma=x)}{dnorm(\hat \mu, \mu = \hat \mu, \sigma=1)}$,

$\frac{|AB|}{|EF|} \cdot dnorm(\hat \mu, \mu = \hat \mu, \sigma=1) = dnorm(A, \mu=\hat \mu, \sigma=x)$

We know everything in the last equation except $x$.

SOLUTION: $\sigma = \sqrt \frac{-(A - \hat \mu) ^ 2} {2 \log \Bigl(\frac{|AB|}{|EF|} \Bigr)}$.

Next question: this estimate is really non-efficient. Can it be improved by choosing different $A$ and averaging?

I attach R code in case if smb wants to play with this estimator. Please, tell me about all improvements you can suggest.

Estimator of SD itself:

find_central_component_mean_and_standard_deviation = function(x) {
  sds <- c()
  density_of_x <- density(x, bw="SJ", kernel="biweight")
  mu = density_of_x$x[which.max(density_of_x$y)]
  EF = max(density_of_x$y)
  density_of_x$y = density_of_x$y / EF * dnorm(0)
  EF = max(density_of_x$y)

  start_to_the_left <- max(which(density_of_x$x < mu - 0.1))
  start_to_the_right <- min(which(density_of_x$x > mu + 0.1))

  previous_value = density_of_x$y[start_to_the_left]
  for (i in seq(from = start_to_the_left, to=1, by=-1)) {
    AB =  density_of_x$y[i]
    if ((AB > previous_value + 10**-10 | AB < 10^-2)) {
    } else {
      if (log(AB / EF) < 0) {
        sds <- c(sds, sqrt(-(density_of_x$x[i] - mu) ** 2 / (2 * log(AB / EF))))
        previous_value = AB

  next_value = density_of_x$y[start_to_the_right]
  for (i in seq(from = start_to_the_right, to=length(density_of_x$x), by=1)) {
    AB =  density_of_x$y[i]
    if ((AB > next_value + 10**-10 | AB < 10^-2)) {
    } else {
      if (log(AB / EF) < 0) {
        sds <- c(sds, sqrt(-(density_of_x$x[i] - mu) ** 2 / (2 * log(AB / EF))))
        next_value = AB
  final_sds <- mean(sds)
  return(c(mu, final_sds))

Simluation procedure:

estimated_sds <- rep(0, 1000)
for (j in 1:1000) {
  vec <- c(rnorm(1000, sd=0.5), runif(1000,min=-6, max=-2), runif(1000,min=2, max=6))
  estimated_sds[j] <- find_central_component_mean_and_standard_deviation(vec)[2]

real_sds <- replicate(1000,sd(rnorm(1000, sd=0.5)))
real_Qns <- replicate(1000,Qn(rnorm(1000, sd=0.5)))

plot(density(na.omit(estimated_sds), bw="SJ"))
lines(density(real_sds, bw="SJ"), col="red")
lines(density(real_Qns, bw="SJ"), col="blue")

LAST UPDATE: https://stat.ethz.ch/pipermail/r-help/2008-October/177121.html - the answer to the problem. Unbiased, accurate, efficient.

  • 4
    $\begingroup$ Because most of the distribution lies beyond the lines you draw, it is not meaningful to declare values beyond those lines to be "outliers"! I'm not questioning your objective, which seems to be to identify some component of a mixture, but am just trying to point out that a different characterization of this objective might be of some assistance to readers in understanding your question and helping you better. In particular, more information about the nature of this distribution and measurement system would be useful. $\endgroup$ – whuber Dec 25 '17 at 23:34
  • 1
    $\begingroup$ Do your "shoulders" contain both "broken" and "unbroken" data points? Or only "broken" ones? And does the "bump" only contain "unbroken" points? If the two distributions do not overlap, you can simply find your two vertical lines, then consider the data between them as a doubly truncated normal distribution and match moments. If you truly have a mixture, you should be able to estimate it, since you already assume the "shoulders" to be uniform. $\endgroup$ – Stephan Kolassa Dec 26 '17 at 9:52
  • $\begingroup$ @whuber I agree, the goal is an extraction of the component, I put the word "outliers" just because my question was marked as "unclear what are you asking" =) I will update the question in 30 mins. $\endgroup$ – German Demidov Dec 26 '17 at 11:41
  • $\begingroup$ @StephanKolassa the bump contain only unbroken points, but the shoulders are mixed. Actually the truncated normal is potentially what I was asking for. I will try to match the moments. Hope it will solve the problem, if not - will update the question. $\endgroup$ – German Demidov Dec 26 '17 at 11:44
  • 1
    $\begingroup$ @StephanKolassa probably I am terribly wrong, but to match the moments of truncated normal with my "bump" I should know much more about the distribution. I need to know $\alpha, \beta, \mu, \sigma$ and can estimate from data only two things (mean and variance) which is not enough to infer these 4 parameters. $\endgroup$ – German Demidov Dec 26 '17 at 12:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.