# Computing standard deviation of discrete normal distribution

I used below pseudocode to generate a discrete normal distribution over 101 points.

mean = 0;
stddev = 1;
lowerLimit = mean - 4*stddev;
upperLimit = mean + 4*stddev;
interval = (upperLimit-lowerLimit)/101;

for ( x = lowerLimit + 0.5*interval ; x < upperLimit; x = x + interval) {
y = exp(-sqr(x)/2)/sqrt(2*PI);
print ("%f %f", x, y);
}
}


When I plot y Vs x I get normal distribution curve as expected. But when I try to calculate standard deviation I use following algorithm (According to http://en.wikipedia.org/wiki/Standard_deviation#Discrete_random_variable)

for i = 1:101
sumsq += y[i]*(x[i]^2)
end
stddev = sqrt(sumsq)


I get $stddev = 3.55$ instead of $1$. Where is the problem?

-
naive question: did you compute the mean as well, to see whether it was indeed 0? (if not, the way you compute your variance is wrong). –  Clement C. Feb 12 '14 at 4:35
@ClementC. yes, I computed mean. It was 2.93714e-15, so nearly zero. –  user13107 Feb 12 '14 at 4:36
Then I reckon it is due to the way you generate the probabilities; in particular, you do not put enoug probability mass near to 0, due to your uniformly spaced grid. One way to fix this would be to get a real probability distribution over the points, by renormalizing all the $y_i$'s (dividing all of them by $\sum_i y_i$) in order to have their sum equal to 1. –  Clement C. Feb 12 '14 at 4:41
@ClementC., this is an algorithmic fix, and the OP would be better off figuring out why it happened that they did not get the total mass equal to 1. There may be weird situations with heavy tailed distributions when any discrete version would have all the moments, yet the continuous version (say Cauchy) won't even have the mean. So one has to be careful with discretized distributions. –  StasK Feb 12 '14 at 5:50

The height of your "density" y should account for the width of the interval that you have discretized your distribution to. In other words, you need to assign the mass $\int_{x_0}^{x_0+h} \phi(x) \, {\rm d}x = \Phi(x_0+h) - \Phi(x_0) \approx \phi(x_0) h$ to the point $x_0$, where $\phi(x)$ is the standard normal density, $\Phi(x)$ is the cdf, and $h$ is your interval. Instead, you assign all of $\phi(x_0)$. You are off by the factor of $h=8/101=0.079$, and when you take the inverse square root of that, you get your 3.55.

-

y isn't a probability distribution.

-
y not? Looks like the usual formula, with exp(-sqr(x)/2) taken as $e^{-x^2/2}.$ –  coffeemath Feb 12 '14 at 4:51
The sum of $y$ is not 1. In fact, for fixed lower and upper limits, increasing the number of points will make sum $y$ go to infinity. You want the area under the normal distribution in intervals centered at the points, and the first and last points being the left and right tails to get a probability distribution. –  Batman Feb 12 '14 at 4:53
I agree, Batman: even though it's the usual formula, the way the x values are generated by the program doesn't make them (approximately) fit a normal distribution. –  coffeemath Feb 12 '14 at 5:08

You need to re-think how you are "discretizing" the normal distribution. You need to either: (1) partition the real line and set the probability of each discrete value to the probabiity of one of these intervals as calculated by the non-discrete version of the normal distribution. or (2) divide the "probabilties" for each value by the sum of the "probability" assigned to all values. This will "renormalize" your values so they sum to 1, and hence represent a probability.

-