How to find the variance(s) of a bivariate normal density such that 95% of the mass is within a certain radius from the mean defined by a point A?

Question

I would like to find the variance of a bivariate normal density (BND), centered at the mean M, such that 95% of its mass is within a certain radius, which depends on the position of a point, A.

(Note: In this toy example, the variances are equal and rho is zero so that the BND will appear as concentric circles.)

The radius is defined as the distance between point M (the mean of the BND, which is fixed) and point A. M is fixed, while A is variable.

In other words, given a certain point A, what is the value of sigma such that 95% of the BND is contained within the radius between M and A?

Since the BND can be decomposed at two univariate densities projecting on two planes, I was wondering if I could just divide the radius by two, and set that as the standard deviation on one dimension. Once I get the standard deviation I can calculate the variances (since they are equal.) In other words, I was wondering if this problem could be the two-dimension version of asking what the variance is if 95% of the density is contained within a certain radius of the mean.

Here is a visualization in r, in which I set the standard deviation as half the radius. It doesn't look right, so I was wondering what would be the correct way to calculate the standard deviation and the variance.

library(mixtools)
library(mnormt)
library(mvtnorm)
library(shape)

# point A
A <- c(3.4, -0.032)
# point M (mean of BND)
M <- c(3.4, 1)
# radius (distance between point and mean of BND)
radius <- M[2] - A[2]

# standard deviation
sd <- radius/2
# sigma
s1 <- sd^2

x.points <- seq(-1,5,length.out=100)
y.points <- seq(-1,5,length.out=100)
z <- matrix(0,nrow=100,ncol=100)

mu1 <- c(3.4,1)
sigma1 <- matrix(c(s1^2,0,0,s1^2),nrow=2)
for (i in 1:100) {
  for (j in 1:100) {
    z[i,j] <- dmvnorm(c(x.points[i],y.points[j]),
                      mean=mu1,sigma=sigma1)
  }
}

par(pty="s")
contour(x.points,y.points,z,xlim=range(2,5), ylim=c(-1,2), nlevels = 5, drawlabels = TRUE)
points(A[1], A[2], pch = 4, col = "red")
plotcircle(mid = c(3.4, 1), r = 0.04, col = "black") 
text(3.65, 0, "A", cex = 1.4)
text(4, 1, "M", cex = 1.4)
text(3.75, 0.4, "radius", cex = 1.2)
segments(3.4, 1, 3.4, 0, col= 'black', cex =1)

^{Created on 2021-03-16 by the reprex package (v0.3.0)}

Answer 1

$$ \sqrt{6}\approx \sqrt{\chi^2_{2, 0.95}} = \sqrt{ (\matrix{r/\sqrt{2} & r/\sqrt{2}})\bigg(\matrix{\sigma^2 & 0\\ 0 &\sigma^2}\bigg)^{-1}\bigg(\matrix{r/\sqrt{2}\\r/\sqrt{2}}\bigg)} $$

This is the Mahalanobis distance in your situation. In more generality, the diagonal elements need not be equal, and the off-diagonal elements are the covariances.

You know your radius $r$. Solve for $\sigma$ or $\sigma^2$.

EDIT

Thinking about this more, I might see some issues if you want correlation or unequal variances, but this does work for your setup and generalizations to higher dimensions.

Answer 2

The OP has revised his question multiple times, and this answer attempts to address an earlier version (highlighted below) which still exists in the current version. My answer below makes no attempt to answer more recent additions to the OP's question which demand to know why the OP's R program outputs (as also provided in the OP's answer to his own question) don't match the (correct) theoretical results provided below and also in the answer by Dave.

Now that the OP's edits have clarified the question somewhat, the question is about a bivariate normal distribution whose joint density ("BND") has circular symmetry about the mean point $(\mu_X,\mu_Y)$. Thus, the contours of the joint density (points ate equal height above the $x$-$y$ plane) are concentric circles whose common center is the mean point $(\mu_X,\mu_Y)$). This tells us that the random variables are independent and have the same variance which is denoted by $\sigma^2$. The joint density has maximum value $\frac{1}{2\pi\sigma^2}$ at $(\mu_X,\mu_Y)$.

So, what should the standard deviation $\sigma$ be so that the probability that $(X,Y)$ lies in the disk of specified radius $r$ centered at $(\mu_X,\mu_Y)$ is $0.95$? Well, if $(X,Y)$ is in the disk, then $(X-\mu_X)^2 + (Y-\mu_Y)^2$ is smaller than $r^2$. But, $W= (X-\mu_X)^2 + (Y-\mu_Y)^2$ is an exponential random variable with rate parameter $\frac{1}{2\sigma^2}$ and so $$P\{(X-\mu_X)^2 + (Y-\mu_Y)^2 > r^2\} = P\{W > r^2\} = \exp\left(-r^2/2\sigma^2\right)$$ which has value $1-0.95 =0.05$ when $$\sigma = \frac{r}{\sqrt{-2 \ln 0.05}} = \frac{r}{\sqrt{2 \ln 20}} = \frac{r}{\sqrt{5.99146\cdots}} = \frac{r}{2.4477\cdots}$$ which is the same result as the one obtained by @Dave but from a purely probabilistic perspective, that is, without getting into Mahalonobis distance and $\chi^2$ random variables and such matters dear to the heart of the statistician.

Turning to the OP's specific question

Since the BND can be decomposed at two univariate densities projecting on two planes, I was wondering if I could just calculate the radius, divide it by two, and set that as the standard deviation on one dimension. Once I get the standard deviation I can calculate the variances (since they are equal.)

Will this guarantee that 95% of the mass is within the circle centered at the mean and with radius equal to the dance (sic) between the mean and A?

it is not clear what is meant by "calculate the radius" since the radius $r$ is given and not something that can be or needs to be calculated. But, dividing the given radius by $2$ and setting it as the standard deviation $\sigma$ makes $\sigma = \frac r2$ (instead of $\frac{r}{2.4477\cdots}$ as described above) and so $$P\{(X-\mu_X)^2 + (Y-\mu_Y)^2 > r^2\} = \exp\left(-r^2/2\sigma^2\right) = \exp(-2) = 0.1353\cdots$$ and so NO, the OP's proposed method does not guarantee that $95\%$ of the probability mass lies in the disk of raduis $r$ centered at $(\mu_X,\mu_Y)$.

The OP has now removed the statement about calculating $r$ which I complained of above. The rest of the latter half of this answer now shows why the OP's revised method does not produce the result the OP hoped for and how much smaller than $r/2$ $\sigma$ must be to achieve the desired $95\%$ probability.

Answer 3

I am having issues implementing the solution proposed by Dilip and Dave.

I agree that for

$$P\{(X-\mu_X)^2 + (Y-\mu_Y)^2 > r^2\} = \exp\left(-r^2/2\sigma^2\right) = 0.05$$

We need $ \sigma=\frac{r}{2.4477\cdots}$

But if I try to implement this in r, and then I compare with the results given by the ellipse() function, I get different results.

I am using the code posted above, with the following changes:

sd <- radius/2.4477

At the end, I add the line:

ellipse(mu=mu1, sigma=sigma1, alpha = .05, npoints = 250, col="red")

As you can see, I get the isoline of the 0.95 quantile of the PDF, but that that not go thorugh point A.

But if I instead divide $\sigma$ by 1.6:

sd <- radius/1.6

And at the end I always run:

ellipse(mu=mu1, sigma=sigma1, alpha = .05, npoints = 250, col="red")