圆形区域上的概率答案

【问题标题】：Probability over a circular region圆形区域上的概率
【发布时间】：2017-10-30 23:31:11
【问题描述】：

假设 z 是一个二元正态（高斯）随机变量，均值为 s 和协方差矩阵 b^2 I_2。我想得到区域 ||x-s|| 上的概率表示某个固定点 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

这是不可能的（因为正方形大于圆形）。我知道有什么问题，但无法查明。你能帮我找出圆形区域上的概率吗？

附： 1) 函数“e2dist”计算两点之间的欧几里得距离。

2) dmvnorm 和 pmvnorm 都来自包“mvtnorm”。

【问题讨论】：

e2dist 和 dmvnorm 的功能来自哪些包？
我相信pmvnorm(lower=c(x[1]-r, x[2]-r), upper=c(x[1]+r, x[2]+r), mean=s, sigma=b*diag(2))[[1]] 返回由lower 和upper 界定的圆形区域中的概率，而不是您问题中建议的方形区域。
@Pascal 我已经编辑了这个问题。抱歉，添麻烦了。函数 e2dist 已在许多空间建模的 r 包中使用。但无论如何我应该给出解释。
@ChrisHolbrook 很抱歉不同意你的观点。 pmvnorm 计算具有任意限制和参数的正态分布函数。浏览页面link。
我知道黎曼积分适用于矩形域。我尝试这种方法只是为了看看它有多接近，因为我不知道计算圆上概率的正确方法。

标签： r probability normal-distribution

【解决方案1】：

您可以使用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 函数返回多变量（不仅是双变量）正态分布的偏移椭圆区域的概率。

【讨论】：

【解决方案2】：

这是一种蛮力方法：

使用问题中的这些：

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

【讨论】：

【解决方案3】：

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

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

【讨论】：