Does standard distance follow the 68-95-99.7 rule?

Question

I'm wanting to do a simple standard distance demonstration for my students in R, but I've come across a conundrum. When I simulate the creation of 10,000 points in a spatial normal distribution, nearly all data fall within 2 standard distances of the mean (far more than the expected ~95%). In fact, it seems as though about 63% falls within 1 standard distance and about 98% falls within 2 standard distances (pretty consistent results through several iterations).

From my understanding, points following a spatial normal distribution also follow the 68-95-99.7 rule with respect to standard distance (source). What's going on?

Here's what I'm using to calculate standard distance and create some visualization:

library(plotrix)

## create function for standard distance
st.dist <- function(x,y) {
  sqrt(sum(((x - mean(x))^2)/length(x)) + (sum((y - mean(y))^2)/length(y)))
}
## assign values to x
x <- rnorm(10000)

## assign values to y
y <- rnorm(10000)

## plot them
plot(x, y,
     xlim = c(mean(x) - 4*st.dist(x,y),
              mean(x) + 4*st.dist(x,y)),
     ylim = c(mean(y) - 4*st.dist(x,y),
              mean(y) + 4*st.dist(x,y)))

## draw standard distances as circles
draw.circle(x = mean(x),
            y = mean(y),
            radius = st.dist(x,y))

draw.circle(x = mean(x),
            y = mean(y),
            radius = 2*st.dist(x,y))


draw.circle(x = mean(x),
            y = mean(y),
            radius = 3*st.dist(x,y))

## create a data frame for the values
df <- data.frame(x=x, y=y, dist=NA)

## calculate distance from mean center to each point
df$dist <- sqrt((mean(x) - df$x)^2 + (mean(y) - df$y)^2)

Checking the percentage of points within 1, 2, and 3 standard distances of the mean center:

> ## percentage of points falling within one standard distance of the mean center
> nrow(df[df$dist<st.dist(x,y),])/nrow(df)
[1] 0.6342
> 
> ## percentage of points falling within two standard distances of the mean center
> nrow(df[df$dist<2*st.dist(x,y),])/nrow(df)
[1] 0.9804
> 
> ## percentage of points falling within three standard distances of the mean center
> nrow(df[df$dist<3*st.dist(x,y),])/nrow(df)
[1] 0.9998
> 

Yet, the x and y values on their own seem to be behaving properly:

## percentage of x coordinates between -1 and 1
> length(x[x<1 & x>-1])/length(x)
[1] 0.6769
> 
> ## percentage of x coordinates between -2 and 2
> length(x[x<2 & x>-2])/length(x)
[1] 0.9545
> 
> ## percentage of x coordinates between -3 and 3
> length(x[x<3 & x>-3])/length(x)
[1] 0.9964
> 
> ## percentage of y coordinates between -1 and 1
> length(y[y<1 & y>-1])/length(y)
[1] 0.6878
> 
> ## percentage of y coordinates between -2 and 2
> length(y[y<2 & y>-2])/length(y)
[1] 0.9514
> 
> ## percentage of y coordinates between -3 and 3
> length(y[y<3 & y>-3])/length(y)
[1] 0.9977
> 

Answer 1

It seems you have defined (from your R code) standard distance (terminology I have not seen before) as the square root of sum of variances along coordinate directions. A more geometrically meaningful definition would be as Mahalanobis distance, but let us stay with your definition. Then

From my understanding, points following a spatial normal distribution also follow the 68-95-99.7 rule with respect to standard distance (source). What's going on?

Short answer: definitely not true. Why? For that we need the distribution of standard distance. So let $(X_1,Y_1), \dotsc, (X_n,Y_n)$ be iid observations from a bivariate normal distribution $\mathcal{bivnorm}(\mu_x,\mu_y,\sigma_x,\sigma_y,\rho)$. Let $D$ be the $n\times 2$ data matrix (with $X_i$, $Y_i$ observations as columns) and $H=I_n - n^{-1} 11^T$ the $n\times n$ centering matrix (a projection matrix), $1=1_n$ is the vector with only $n$ 1's. Then we can write $\text{stdist}= \sqrt{\text{trace}(n^{-1} D^T H D))}$. So what is the distribution? I will come back for a full (I hope) answer, but for now some references: Traces and cumulants of quadratic forms in normal variables.

But for the very specific case $\rho=0$ this is simple, then we have $$ \text{stdist}= \sqrt{\frac{n-1}{n}(S_x^2 + S_y^2)}=\sqrt{\frac{1}{n}}\sqrt{(n-1)S_x^2 + (n-1)S_y^2} $$ and the sample variances are in this case independent, so the sum under the square root sign in the last factor have a chisquared distribution with $2n-2$ degrees of freedom. That will answer your question in that case.