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'.
You can get this probability with the shotGroups package:
> library(shotGroups)
> pmvnEll(r=0.5, sigma=diag(2), mu=c(0.1,0.1), e=diag(2), x0=c(0,0))
[1] 0.1164051
More generally, the pmvnEll function returns the probability of an offset ellipsoidal region for a multivariate (not only bivariate) normal distribution.