How to calculate 2D standard deviation, with 0 mean, bounded by limits

Question

My problem is as follows: I drop 40 balls at once from a certain point, a few meters over the floor. The balls roll, and comes to a rest. Using computer vision, I calculate the center of mass in the X-Y plane. I am only interested in the distance from the center of mass to each ball, which is calculated using simple geometry. Now, I want to know the one-sided standard deviation from the center. So, I would be able to know that a certain number of balls are within one std radius, more balls within 2*std radius and so on. How do I calculate the one-sided standard deviation? A normal approach would state that half of the balls be on the "negative side" of 0 mean. This of course makes no sense in this experiment. Do I have to make sure that the balls conform to the standard distribution? Thank you for any help.

Answer 1

To characterize the amount of 2D dispersion around the centroid, you just want the (root) mean squared distance,

$$\hat\sigma=\text{RMS} = \sqrt{\frac{1}{n}\sum_i\left((x_i - \bar{x})^2 + (y_i - \bar{y})^2\right)}.$$

In this formula, $(x_i, y_i), i=1, 2, \ldots, n$ are the point coordinates and their centroid (point of averages) is $(\bar{x}, \bar{y}).$

---

The question asks for the distribution of the distances. When the balls have an isotropic bivariate Normal distribution around their centroid--which is a standard and physically reasonable assumption--the squared distance is proportional to a chi-squared distribution with two degrees of freedom (one for each coordinate). This is a direct consequence of one definition of the chi-squared distribution as a sum of squares of independent standard normal variables, because $$x_i - \bar{x} = \frac{n-1}{n}x_i - \sum_{j\ne i}\frac{1}{n}x_j$$ is a linear combination of independent normal variates with expectation $$\mathbb{E}[x_i - \bar{x}] = \frac{n-1}{n}\mathbb{E}[x_i] -\sum_{j\ne i}\frac{1}{n}\mathbb{E}[x_j] = 0.$$ Writing the common variance of the $x_i$ as $\sigma^2$, $$\mathbb{E}[\left(x_i -\bar{x}\right)^2]=\text{Var}(x_i - \bar{x}) = \left(\frac{n-1}{n}\right)^2\text{Var}(x_i) + \sum_{j\ne i}\left(\frac{1}{n}\right)^2\text{Var}(x_j) = \frac{n-1}{n}\sigma^2.$$ The assumption of anisotropy is that the $y_j$ have the same distribution as the $x_i$ and are independent of them, so an identical result holds for the distribution of $(y_j - \bar{y})^2$. This establishes the constant of proportionality: the squares of the distances have a chi-squared distribution with two degrees of freedom, scaled by $\frac{n-1}{n}\sigma^2$.

The most severe test of these equations is the case $n=2$, for then the fraction $\frac{n-1}{n}$ differs the most from $1$. By simulating the experiment, both for $n=2$ and $n=40$, and overplotting the histograms of squared distances with the scaled chi-squared distributions (in red), we can verify this theory.

Each row shows the same data: on the left the x-axis is logarithmic; on the right it shows the actual squared distance. The true value of $\sigma$ for these simulations was set to $1$.

These results are for 100,000 iterations with $n=2$ and 50,000 iterations with $n=40$. The agreements between the histograms and chi-squared densities are excellent.

---

Although $\sigma^2$ is unknown, it can be estimated in various ways. For instance, the mean squared distance should be $\frac{n-1}{n}\sigma^2$ times the mean of $\chi^2_2$, which is $2$. With $n=40$, for example, estimate $\sigma^2$ as $\frac{40}{39}/2$ times the mean squared distance. Thus an estimate of $\sigma$ would be $\sqrt{40/78}$ times the RMS distance. Using values of the $\chi^2_2$ distribution we can then say that:

• Approximately 39% of the distances will be less than $\sqrt{39/40}\hat\sigma$, because 39% of a $\chi^2_2$ distribution is less than $1$.

• Approximately 78% of the distances will be less than $\sqrt{3}$ times $\sqrt{39/40}\hat\sigma$, because 78% of a $\chi^2_2$ distribution is less than $3$.

And so on, for any multiple you care to use in place of $1$ or $3$. As a check, in the simulations for $n=40$ plotted previously, the actual proportions of squared distances less than $1, 2, \ldots, 10$ times $\frac{n-1}{n}\hat\sigma^2$ were

0.3932 0.6320 0.7767 0.8647 0.9178 0.9504 0.9700 0.9818 0.9890 0.9933

The theoretical proportions are

0.3935 0.6321 0.7769 0.8647 0.9179 0.9502 0.9698 0.9817 0.9889 0.9933

The agreement is excellent.

---

Here is R code to conduct and analyze the simulations.

f <- function(n, n.iter, x.min=0, x.max=Inf, plot=TRUE) {
  #
  # Generate `n.iter` experiments in which `n` locations are generated using
  # standard normal variates for their coordinates.
  #
  xy <- array(rnorm(n*2*n.iter), c(n.iter,2,n))
  #
  # Compute the squared distances to the centers for each experiment.
  #
  xy.center <- apply(xy, c(1,2), mean)
  xy.distances2 <- apply(xy-array(xy.center, c(n.iter,2,n)), c(1,3), 
                         function(z) sum(z^2))
  #
  # Optionally plot histograms.
  #
  if(plot) {
    xy.plot <- xy.distances2[xy.distances2 >= x.min & xy.distances2 <= x.max]

    hist(log(xy.plot), prob=TRUE, breaks=30,
         main=paste("Histogram of log squared distance, n=", n),
         xlab="Log squared distance")
    curve(dchisq(n/(n-1) * exp(x), df=2) * exp(x) * n/(n-1), 
          from=log(min(xy.plot)), to=log(max(xy.plot)), 
          n=513, add=TRUE, col="Red", lwd=2)

    hist(xy.plot, prob=TRUE, breaks=30,
         main=paste("Histogram of squared distance, n=", n),
         xlab="Squared distance")
    curve(n/(n-1) * dchisq(n/(n-1) * x, df=2), 
          from=min(xy.plot), to=max(xy.plot), 
          n=513, add=TRUE, col="Red", lwd=2)  
  }
  return(xy.distances2)
}
#
# Plot the histograms and compare to scaled chi-squared distributions.
#
par(mfrow=c(2,2))
set.seed(17)
xy.distances2 <- f(2, 10^5, exp(-6), 6)
xy.distances2 <- f(n <- 40, n.iter <- 50000, exp(-6), 12)
#
# Compare the last simulation to cumulative chi-squared distributions.
#
sigma.hat <- sqrt((n / (2*(n-1)) * mean(xy.distances2)))
print(cumsum(tabulate(cut(xy.distances2, 
                    (0:10) * (n-1)/n * sigma.hat^2))) / (n*n.iter), digits=4)
print(pchisq(1:10, df=2), digits=4)

Answer 2

I think you have some things a bit confused. It's true that distance can't be negative, but that doesn't affect calculation of the standard deviation. Although it means the distribution of distances can't be exactly normal, it could still be close; but even if its is far from normal, there is still a standard deviation.

Also, there is no "one sided" standard deviation - you may be thinking of hypothesis tests (which can be one sided or two sided). In your title, you say mean is 0, but the mean distance won't be 0 (unless the balls are in a stack 40 balls high!) and you say there are limits - there could be limits, if the balls are dropped in a room then they can't be farther from the center than the distance to the nearest wall. But unless some of the balls bounce against a wall, that won't affect things.

So, once you have the 40 distances you calculate the standard deviation (and mean, median, interquartile range, etc) using standard methods. You can also make plots of the distance (e.g. quantile normal plot, box plot) to see if it is roughly normally distributed (if that's of interest).

Answer 3

Its been a while since this was asked, but the answer to the question is that this is the 2D distribution named the Rayleigh distribution. Here the assumption is that the Rayleigh shape factor is equal to both the standard deviations of the X and Y coordinates. In practice the value of the shape factor would be calculated from the pooled average of the standard deviation of X and Y.

starting with $$ X \sim \mathcal{N}(\mu_x,\sigma_x^2)$$, and $$Y \sim \mathcal{N}(\mu_y,\sigma_y^2)$$

use bivariant normal distribution. $$ f(x,y) = \frac{1}{2 \pi \sigma_x \sigma_y \sqrt{1-\rho^2}} \exp\left( -\frac{1}{2(1-\rho^2)}\left[ \frac{(x-\mu_x)^2}{\sigma_x^2} + \frac{(y-\mu_y)^2}{\sigma_y^2} - \frac{2\rho(x-\mu_x)(y-\mu_y)}{\sigma_x \sigma_y} \right] \right)$$

translate to point $$(\mu_x, \mu_y)$$ and assume $$\rho = 0$$.

Also assume that $$\sigma_x^2 = \sigma_y^2$$ so replace both with $$\sigma^2$$

then the 2-D distribution is expressed as the radius around point $$(\mu_x, \mu_y)$$ which is known as the Rayleigh distribution.

$$PDF(r; \sigma) = \frac{r}{\sigma^2 } \exp\left( - \frac{r^2}{2\sigma^2} \right) $$ where $$\sigma = \sigma_x = \sigma_y$$ and $$r_i = \sqrt{(x_i - \mu_x)^2 + (y_i - \mu_y)^2}$$

$$ CDF(r; \sigma) = 1 - \exp\left( - \frac{r^2}{2\sigma^2} \right)$$

Of course this is for the continuous distribution. For a sample of just 40 balls there is no exact solution. You'd need to do a Monte Carlo Analysis with a sample of 40 balls. Taylor, M. S. & Grubbs, Frank E. (1975). "Approximate Probability Distributions for the Extreme Spread" found estimates for the Chi distribution and the log-normal for that would fit the distribution of a sample.

---

Edit - Despite Wuber's doubt, the theoretical proportions he calculated are:

0.3935 0.6321 0.7769 0.8647 0.9179 0.9502 0.9698 0.9817 0.9889 0.9933

From the CDF function the cumulative Sigma values for r (in sigmas) equal to range from:

0-1, 0-2, 0-3, ... , 0-10

are:

0.3935, 0.6321, 0.7769, 0.8647, 0.9179, 0.9502, 0.9698, 0.9817, 0.9889, 0.9933

Answer 4

The normal distribution, both positive and negative values, makes sense if you recognize that this normal distribution is for radius or "distance from centroid". The other variable, angle, is random and is uniformly distributed from 0-pi