# How to estimate mean and standard deviation of a normal distribution from noisy data?

I have $n$ observations, $x_i$ following a normal distribution. I would like to estimate $\mu$ and $\sigma$ from my samples. Normally I would simply estimate

$\mu=(\sum x_i)/n$ and
$\sigma^2=\sum (x_i-\mu)^2/(n-1)$.

However I know that my samples $x_i$ are noisy. I know that the 95% confidence interval of

$x_i$ is $[x_i-t, x_i+t]$

with the same $t$ for all $x_i$. How do I take the noise into account when estimating $\sigma$?

If I don’t take the confidence interval into account the estimation of $\sigma$ might be too large.

• What is a confidence interval for $x_i$? The $x_i$ are data, immutable and unchangeable as far as you are concerned; you can have 100% confidence that the value of $x_i$ is exactly $x_i$, no more, no less. Jan 22, 2015 at 12:27
• You are right. Thanks for helping me to state my question more precisely. The x_i are estimations an unknown value X_i. I can measure X_i with a certain device. The device makes a certain error. I assume that the measurements are normally distributed. I have (unfortunately) only two measurements v_1 and v_2 with the device to determine X_i. So I estimate X_i\app x_i=(v_1+v_2)/2. So I know that X_i\in [x_i-t, x_i+t] with a probability of 95% with the usual t for confidence intervals. Since I have only 2 measurements my confidence interval will be large. Jan 22, 2015 at 13:04
• Please edit the information into your question. Jan 22, 2015 at 14:02
• Is this not the prototypical ANOVA (components of variance) problem? Do you have two replicate measurement for each $x$ ("balanced" data) or have you replicated only one of them ("unbalanced" data)?
– whuber
Jan 22, 2015 at 17:30

Partially answered in comments: Is this not the prototypical ANOVA (components of variance) problem? Do you have two replicate measurement for each x ("balanced" data) or have you replicated only one of them ("unbalanced" data)? – whuber

Since the OP didn't respond to the implicit questions by whuber, but I can show (with simulated data) how this could be done in R when there are replicate measurements for each $$X_i$$:

library(lme4)

set.seed(7*11*13) # my public seed
X <- 1:10
r <- 3
mu <- rep(X,rep(3,length(X)))
obs <- rnorm(length(mu),mu,1)

ind <- as.factor(mu)
mydata <- data.frame(obs, mu, ind)
mod <- lmer(obs ~ 1 + ind + (1 | ind), data=mydata )

mod
Linear mixed model fit by REML ['lmerMod']
Formula: obs ~ 1 + ind + (1 | ind)
Data: mydata
REML criterion at convergence: 74.71
Random effects:
Groups   Name        Std.Dev.
ind      (Intercept) 3.14
Residual             1.19
Number of obs: 30, groups:  ind, 10
Fixed Effects:
(Intercept)         ind2         ind3         ind4         ind5         ind6
1.6086      -0.6777       1.2451       2.5058       3.8129       4.7079
ind7         ind8         ind9        ind10
5.1083       7.2207       6.7395       7.7469


If you really know the true values for your specimens (as measured by some much more precise instrument, say) you can do:

 mod2 <- lmer(obs ~ 0 + offset(mu) + (1 | ind), data=mydata )
mod2
Linear mixed model fit by maximum likelihood  ['lmerMod']
Formula: obs ~ 0 + offset(mu) + (1 | ind)
Data: mydata
AIC      BIC   logLik deviance df.resid
96.7257  99.5281 -46.3629  92.7257       28
Random effects:
Groups   Name        Std.Dev.
ind      (Intercept) 0.000
Residual             1.135
Number of obs: 30, groups:  ind, 10
No fixed effect coefficients