如何从超椭球表面均匀采样(恒定马氏距离)?

机器算法验证 随机生成 均匀分布
2022-03-12 09:29:33

在实值多元情况下,有没有办法从表面上均匀采样点,其中马氏距离与平均值的距离是常数?

编辑:这只是归结为从满足方程的超椭球表面均匀采样点,

(xμ)TΣ1(xμ)=d2.

更准确地说,“均匀地”是指采样使得每个区域元素dA的超曲面包含相同的概率质量。

1个回答

当不同的椭球轴没有太大差异时,使用拒绝采样是可行的(差异很大,您拒绝很多,使其不太可行)

  • (1) 超球面上的样本
  • (2) 将其挤压成超椭球体
  • (3) 计算表面积被挤压的速率
  • (4) 按该比率拒绝样品。

二维示例

例子

set.seed(1)
#some matrix to transform n-sphere (in this case 2x2)
m <- matrix(c(1, 0.55, 0.55, 0.55), 2)

# sample multinomial with identity covariance matrix
x <- cbind(rnorm(3000, 0, 1), rnorm(3000, 0, 1))
l1 <- sqrt(x[,1]^2 + x[,2]^2)

# perpendicular vector
per <- cbind(x[,2], -x[,1])

# transform x
x <- x %*% m
# transform perpendicular vector (to see how the area transforms)
per2 <- per %*% m

# get onto unit-"sphere"/ellipsoid
x <- x/l1

# this is how the area contracted
contract <- sqrt(per2[,1]^2 + per2[,2]^2) / sqrt(per[,1]^2 + per[,2]^2)

# then this is how we should choose to reject samples 
p <- contract/max(contract)

# rejecting
choose <- which( rbinom(n=length(p), size=1, p=p) == 1)

#plotting
plot(x[1:length(choose), 1], x[1:length(choose), 2],
     xlim=c(-1.2, 1.2), ylim=c(-1.2, 1.2),
     xlab = expression(x[1]), ylab = expression(x[2]),
     bg=rgb(0, 0, 0, 0.01), cex=0.6, pch=21, col=rgb(0, 0, 0, 0.01))
title("squeezed uniform circle \n ")

#plotting
plot(x[choose,1], x[choose,2],
     xlim=c(-1.2, 1.2), ylim=c(-1.2, 1.2),
     xlab = expression(x[1]), ylab = expression(x[2]),
     bg=rgb(0, 0, 0, 0.01), cex=0.6, pch=21, col=rgb(0, 0, 0, 0.01))
title("squeezed uniform circle \n  with rejection sampling")