离散逼近二元正态分布

Question

我希望进行二元正态分布的离散逼近。也就是说，我希望计算一个矩阵，其中每个条目都是掉落到下图中小方块内的概率。

以下是我目前做的内容：

library(mvtnorm)
library(graphics)

euclide = function(x,y){sqrt(x^2+y^2)}
maxdist = 40
sigma = diag(2)
m = matrix(0,ncol=maxdist*2 + 1, nrow=maxdist*2 + 1)
for (row in -maxdist:maxdist){
    for (col in -maxdist:maxdist){
        if ( euclide(abs(row), abs(col)) < maxdist ){
            lower = c(row-0.5, col-0.5)
            upper = c(row+0.5, col+0.5)
            p = pmvnorm(lower = lower , upper = upper, mean = c(0,0), sigma = sigma)    
        } else {
            p = 0
        }
        m[row + maxdist + 1,col + maxdist + 1] = p
    }
}
m = m[rowSums(m)!=0,colSums(m)!=0]
contour(m, levels = exp(-20:0), xlim=c(0.3,0.7), ylim=c(0.3,0.7))

这个程序可以正常工作。但是当maxdist较大时，速度会比较慢。我希望能够提高计算速度，但这不是我的主要问题...

我的主要问题是，使用我的方法时，我无法改变靠近中心的小方块数量，以便更好地逼近平均值。我只能在周围添加方块。换句话说，我想能够设置双变量正态分布的两个轴的方差。

- Remi.b

---

你听说过Tauchen吗？以及如何从N(0,1)的随机数生成联合正态变量？（例如，这里的幻灯片269-272页）或者你是为了其他特定目的而进行这个练习吗？ - MichaelChirico

请问您构建矩阵的目的是什么？ - Robert Dodier

这是一个分散核（到给定距离的单元的分散概率），由于单元内有子单元，因此细节略微复杂。模拟存在于单元中并随机交配的个体比存在于空间连续体中更快计算。这就是为什么我想要这种离散逼近的原因。 - Remi.b

@Remi.b 好的，谢谢你提供的信息。我在想这个关键操作是否会是卷积（将该转移矩阵与当前状态分布进行卷积，以得到下一个状态的分布）。如果是这样，请考虑通过快速傅里叶变换来计算卷积，以提高效率。 - Robert Dodier

Answer 1

我不是R语言专家，但我确定正态分布有一个CDF函数。如果您想要的确实是每个方块落入的概率矩阵，我们可以使用这个CDF函数来得到答案。由于二维正态分布具有独立的边缘分布，因此这里的问题只是针对由轴位置[x_left，x_right]和[y_left，y_right]描述的每个方块提出两个问题：

1. 一维正态随机变量在区间[x_left，x_right]内的概率是多少？ 2. 另一个独立的一维正态随机变量在区间[y_left，y_right]内的概率是多少？

由于这两个变量是独立的，因此该方块的完整概率为：

P = (CDF(x_right) - CDF(x_left))*(CDF(y_right) - CDF(y_left))

这是一个精确的答案，计算时间不应该是一个问题！编辑：我还应该说，您可以选择一个在零附近每个轴上有更多刻度的网格，以获得所需的分辨率。上面每个正方形的概率公式仍然有效。

- Daniel Johnson

Answer 2

下面是一个简单的实现。就像@DanielJohnson所说的，您可以使用单变量正态分布的cdf形式，但它应该与使用pmvnorm相同，如下所示。使用pnorm的版本要快得多。

## Choose the matrix dimensions
yticks <- xticks <- seq(-3, 3, length=100)
side <- diff(yticks[1:2])  # side length of squares
sigma <- diag(2)               # standard devs. for f2
mu <- c(0,0)                # means

## Using pnorm
f <- Vectorize(function(x, y, side, mu1, mu2, s1, s2)
    diff(pnorm(x+c(-1,1)*side/2, mu1, s1)) * diff(pnorm(y+c(-1,1)*side/2, mu2, s2)),
    vec=c("x", "y"))

## Using pmvnorm
f2 <- Vectorize(function(x, y, side, mu, sigma)
    pmvnorm(lower=c(x,y)-side/2, upper=c(x,y)+side/2, mean=mu, sigma=sigma),
                vec=c("x", "y"))

## get prob. of squares, mu are means, s are standards devs.
mat <- outer(xticks, yticks, f, side=side, mu1=0, mu2=0, s1=1,s2=1)
mat2 <- outer(xticks, yticks, f2, side=side, mu=mu, sigma=sigma)

## test equality
all(abs(mat2-mat) < 1e-11)  # TRUE
all.equal(mat2, mat)        # TRUE

## See how it looks
library(lattice)
persp(mat, col="lightblue", theta=35, phi=35, shade=0.1)

- Rorschach