2
$\begingroup$

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)

$\endgroup$
10
  • 1
    $\begingroup$ It will have to do with something call Mahalanobis distance. $\endgroup$
    – Dave
    Mar 9, 2021 at 17:40
  • 1
    $\begingroup$ The radius (divided by the standard deviation of either coordinate) has a $\chi(2)$ distribution, so just pick its 95th percentile, which is close to $\sqrt{6}.$ $\endgroup$
    – whuber
    Mar 9, 2021 at 17:47
  • $\begingroup$ @whuber, if I understand correctly, if the ratio between the radius and the standard deviation is approximately equal to 2.44, then, to get the standard deviation, should I divide the radius by 2.44 instead of 2 as I was thinking initially? $\endgroup$
    – Emy
    Mar 9, 2021 at 17:58
  • $\begingroup$ @whuber The claim that a circle of (standardized) radius approximately $\sqrt{6}$ about an arbitrary point in the plane contains $95\%$ of the probability mass from a standard bivariate normal distribution of independent random variables is not correct. There seen to be multiple interpretations of what the OP is asking. $\endgroup$ Mar 10, 2021 at 15:13
  • $\begingroup$ @Dilip I did not make such a claim and find it difficult to interpret anything in the question as suggesting that. I will agree that I might lack the imagination to conceive of that interpretation, but the repeated references to a "radius" around a "point" seem pretty clear to me. Are you taking the density to be fixed, the point $A$ to be arbitrary, and then asking for the minimal radius of a disk centered at $A$ containing 95% of the probability? $\endgroup$
    – whuber
    Mar 10, 2021 at 15:37

3 Answers 3

3
$\begingroup$

$$ \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.

$\endgroup$
3
  • $\begingroup$ To make sure you understand what's going on with Mahalanobis distance here, write out the equation for your same situation but with three independent variables. $\endgroup$
    – Dave
    Mar 9, 2021 at 18:36
  • $\begingroup$ thanks a lot, by the third independent variable you mean \rho, right? $\endgroup$
    – Emy
    Mar 9, 2021 at 18:38
  • 1
    $\begingroup$ @Emy I mean a trivariate distribution with three independent Gaussian variables. Do it in the full generality with the covariances if you want, but that makes the algebra messy without doing much to show that you understand what's going on with Mahalanobis distance. $\endgroup$
    – Dave
    Mar 9, 2021 at 18:40
3
$\begingroup$

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.

$\endgroup$
10
  • $\begingroup$ I get $r/\sigma = \sqrt{\chi^2_{2, 0.95}}\approx\sqrt{5.99}\approx\sqrt{2\ln 20}$, same as you. $\endgroup$
    – Dave
    Mar 10, 2021 at 11:31
  • $\begingroup$ @Dilip Thank you for the answer. The only correction I would make is that, in my description, the joint density is not centered at A. Point A is the point that defines the radius as a distance from the joint density mean M. Apart from this, I have a question: does this mean that to find the standard deviation, given a radius, I need to divide the radius by the square root of 6? That would mean dividing by 2.44 instead of 2 as I had thought initially. $\endgroup$
    – Emy
    Mar 10, 2021 at 12:00
  • $\begingroup$ Large numbers of samples? What I did is correct at the population level without doing any sampling at all. $\endgroup$
    – Dave
    Mar 10, 2021 at 15:16
  • $\begingroup$ @Dilip, in the question, I have edited "calculate the radius" to "divide the radius" since as you point out, once we know where point A is, the radius is fixed. $\endgroup$
    – Emy
    Mar 11, 2021 at 15:11
  • 1
    $\begingroup$ @Emy I answered the question you asked. If you now have a different question, please post it is a new question. If someone comes along with the question you had, I want them to see my (correct) answer. Likewise, if someone sees my answer to you updated question, I do not want them to be misled by my answer to a different question than the one you originally asked. Please do link this question and say that it’s a follow-up, but moving the goalposts after you’ve gotten answers is discouraged. $\endgroup$
    – Dave
    Mar 16, 2021 at 12:04
-2
$\begingroup$

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")

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.