I have been bugged by a simple problem in statistics recently. Let's assume that I have made a set of measurements of a certain quantity, each with an uncertainty estimate. I have a set of tuples $\{x_i,s_i\}$ ($x$ is a value, $s$ its uncertainty).
The goal is to calculate the most likely true value $y$ and its uncertainty $S$.
My solution: I assume that each measurement $x_i$ follows a Gaussian distribution with a standard deviation $s_i$, meaning that the probability that $y$ is the true value provided $x_i$ is $$ P(y|x_i,s_i) \propto \exp \left ( -\frac{(x_i-y)^2}{2s_i^2} \right ) \,. $$ Therefore the total probability (the measurements are independent) that $y$ is the true value should be $$ P(y|\{x_i,s_i\}) \propto \exp \left ( -\sum_i \frac{(x_i-y)^2}{2s_i^2} \right ) \,. $$ Searching for the most probable $y$, I just set $\partial P / \partial y = 0$. I then get $$ \sum_i \frac{y}{s_i^2} = \sum_i \frac{x_i}{s_i^2} \, . $$ Therefore the most probable value should just be an average measurements, using the squares of reciprocals of their uncertainties as weights. So far so good.
I would then like to know that the standard deviation (its estimate) in $y$ is. I therefore look at $P(y+a|\{x_i,s_i\})/P(y|\{x_i,s_i\})$ (for some $a$).
What I get is $$ \frac{P(y+a)}{P(y)} = \exp \left ( - \sum_i \frac{a^2}{2s_i^2} \right ) \, , $$ therefore the standard deviation should be $$ \frac{1}{S^2} = \sum_i \frac{1}{s_i^2} \, . $$
That all makes sense. Provided that all the errors are the same, then $S = s/\sqrt{N}$. However, if I set all errors to 0 ($s_i = 0$), then the $S$ is also be 0. But there is still some uncertainty due to the spread of values $\{x_i\}$.
My question is: How should I approach the problem in order to incorporate both the spread of $\{x_i\}$ and their uncertainties $\{s_i\}$ in the total uncertainty $S$?
Thanks a lot.
SSF
The way this usually goes is slightly different from what you've written. Usually we assume we have a bunch of independent Gaussian variables $X_1,\dots,X_n$ all of which have mean $\mu$ and which have (perhaps unknown) standard deviations $\sigma_i$. If we know $\sigma_i$ then the least squares estimator of $\mu$ is
$$\sum_{k=1}^n \frac{X_i}{\sigma_i^2}.$$
It is central here that the means are assumed to be $\mu$. In scientific applications we would say that $X_i$ have only random error. If their means deviate from $\mu$ at all, that deviation is systematic error, which pure probability modeling is not very good at handling.
If you find that your actual sample numbers $x_i$ are highly spread out relative to the values of $\sigma_i$ that you have recorded (for instance, if all of the $\sigma_i$ are $1$ and your values range from $-10000$ to $10000$), then you have reason to believe that your assumed values of $\sigma_i$ are incorrect. But this is an issue at the level of modeling, not at the level of mathematics. The mathematics just tells you that if your model is correct, then such a spread is extremely unlikely.