圆形区域的概率

Question

假设 z 是双变量正态（高斯）随机变量，均值为 s，协方差矩阵为 b^ 2 I_2。我想获得 ||x-s|| 区域的概率<= r 对于一些不动点 x 和常量 r > 0。我正在使用 R 软件来计算它。下面是我尝试估计概率的方法的例子 -

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

当我计算正方形区域 (x-c(r,r)) 到 (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

这是不可能的（因为正方形比圆大）。我知道有问题，但找不到。你能帮我求出圆形区域的概率吗？

P.S。 1) 函数 'e2dist' 计算两点之间的欧式距离。

2) dmvnorm 和 pmvnorm 都来自软件包 'mvtnorm'。

Answer 1

这是一种蛮力方法：

使用问题中的这些：

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)

从多元正态分布中抽取大量样本并计算这些样本在感兴趣区域中的比例

y <- rmvnorm(1000000,mean=s, sigma=b*diag(2))

#proportion of mvn distn in circular region (radius r) centered at x
dyx <- e2dist(x,y) #distances between y and x
mean(dyx < r)
>[1] 0.117238

#proportion of mvn distn in circular region (radius r) centered at s
dys <- e2dist(s,y) #distances between y and s
mean(dys < r)
>[1] 0.118308

对于方形区域，结果与pmvnorm非常一致，但可能需要非常大量的随机样本

#proportion of mvn distn in square region
mean( y[,1] >= -r & y[,1] <= r &
      y[,2] >= -r & y[,2] <= r)
>[1] 0.145965

#compare to...
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

Answer 2

你可以用shotGroups包得到这个概率：

> 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

更一般地说，pmvnEll 函数 returns 多变量（不仅是双变量）正态分布的偏移椭圆区域的概率。

Answer 3

扩展我的评论：您的方法原则上是正确的，修正两个错误会产生近似正确的结果：

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(x,grid.pts)<=r                               # circle is centered around x not s
tt=which(ttt==T, arr.ind=T)
circle.in.pts=grid.pts[tt,]
mean(dmvnorm(circle.in.pts,s,b*diag(2))) * pi*r^2         # need to multiply by area

# Output:
# [1] 0.1164057

圆形区域的概率

Probability over a circular region

r

normal-distribution

probability