Variance of noisy data makes unexpected jump I observer a somewhat surprising result.
I have a probability distribution of some object in space. I would like to find the variance of that probability distribution. However, there is a catch. The probability distribution I get is corrupted by small amounts of white noise. My naive approach to fix it is as follows:


*

*Subtract minimum from the corrupted function

*Normalize it

*Calculate the variance


Somewhat surprisingly, for very mild amounts of noise my estimator can be wrong by a factor of 10 !!!
Why so much, and is there a better way to estimate variance in my case? For my actual application I do not know the variance of the noise.
def var(x, y):
    p = y - np.min(y)
    p /= np.sum(y)
    mu = p.dot(x)
    return p.dot((x - mu)**2)


gau = lambda mu, s2: np.exp(-(mu**2)/s2/2) / np.sqrt(2 * np.pi * s2)

x = np.linspace(0, 100, 100)
y = gau(x - 50, 50)
ynoisy = y + np.random.normal(0, 0.001, 100)

print("Original variance", var(x, y))
print("Noisy variance", var(x, ynoisy))

plt.figure()
plt.plot(x, y)
plt.plot(x, ynoisy)
plt.show()


Original variance 49.99999993042545
Noisy variance 412.4350189819331


 A: Ok, now it is much clearer what you want to achieve and also why the proposed approach doesn't work. This is because the minimum of the noisy function can occur anywhere, most likely in the outer flanks where the function is close to zero. So imagine a large negative minimum to occur close to the ends of the plot. This means that you will create a huge offset in order to produce a pdf for which to compute the variance. Unlike a gaussian, this offset (together with the remoteness of the minimum from the maximum) then translates into a large probability mass in those outer regions, which will heavily distort your variance measure. 
So how can that be remedied? I suggest you to try what is called Laplaces method (https://en.wikipedia.org/wiki/Laplace%27s_method), which boils down to fitting a second-order polynomial to the logarithm of your noisy function (or equivalently to the local approximation of the function by a gaussian). 
For this approach to work however you need to make sure that (a) the function can reasonably be approximated locally by a gaussian 
and 
(b) that this local region consists of positive values only
