Suppose *z* is a bivariate normal (Gaussian) random variable with mean at *s* and covariance matrix *b^2 I\_2*. I want to get the probability over the area *||x-s|| <= r* for some fixed point *x* and constant *r > 0*. I was using R software to compute it. The following is an example of the way I tried to estimate the probability -

```
library(mvtnorm)
e2dist=function(x,y) # x is vector, y is matrix
{
  a=sqrt((x[1]-y[,1])^2 + (x[2]-y[,2])^2)
  return(a)
}
r=0.5
b=1
s=c(0.1,0.1)
x=c(0,0)
a1=seq(x[1]-r, x[1]+r, length.out=1000)
a2=seq(x[2]-r, x[2]+r, length.out=1000)
grid.pts=as.matrix(expand.grid(a1,a2))
ttt=e2dist(s,grid.pts)<=r
tt=which(ttt==T, arr.ind=T)
circle.in.pts=grid.pts[tt,]
mean(dmvnorm(circle.in.pts,s,b*diag(2)))
> [1] 0.1503632

```

This estimate of the probability is not right as when I compute the true probability over the square area (x-c(r,r)) to (x+c(r,r)) 

```
pmvnorm(lower=c(x[1]-r, x[2]-r), upper=c(x[1]+r, x[2]+r), mean=s, sigma=b*diag(2))[[1]]
> [1] 0.1452895

```

which is not possible (as the square is larger than the circle). I know there is something wrong, but could not find out. Can you help me to find the probability over the circular area ?

P.S. 1) The function 'e2dist' computes the euclidean distance between two points. 

2) Both dmvnorm and pmvnorm are from package 'mvtnorm'.