2015-07-04 128 views
3

我想提出的离散近似的二元正态分布。也就是说,我想计算一个矩阵,其中每个条目都是落入下图中一个小平方的概率。离散逼近二元正态分布

enter image description here

这里是我做的那么远。

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)) 

enter image description here

它工作正常。这很慢(对于大的maxdist),但我希望能够提高计算时间。但是,这是不是我的主要问题...

的主要问题是,我的方法,我不能改变接近中心做出更好的近似接近均值的小方块的数量。我只能在周围添加方块。换句话说,我希望能够设置双变量正态分布的两个轴的方差。

+0

你听说过[Tauchen]的(http://www.sciencedirect.com/science/article/pii/0165176586901680)?从'N(0,1)'生成联合正态变量画? (例如[slide 269-272 here](http://www.ssc.upenn.edu/~fdiebold/Teaching706/TimeSeriesSlides.pdf))或者你是否正在为另一个特定目的而进行此练习? – MichaelChirico

+0

我可以问你想用什么目的来使用你正在构建的矩阵? –

+0

这是一个分散的核心(分散到一个给定距离的细胞的概率),细节的细节稍微复杂一些,因为细胞内有子细胞。模拟在细胞中存在的个体和随机交配比在空间连续体上存在的个体在计算上更快。这就是为什么我想要这种离散近似。 –

回答

2

下面是一个简单的实现。像@丹尼尔约翰逊说,你可以使用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) 

enter image description here

3

我不是的R的人,但我敢肯定有正态分布一个CDF功能。如果你想要的是从字面上看落入每个方格的概率矩阵,我们可以使用这个CDF函数来得到答案。由于2D正态分布具有独立的边缘分布,问题这里只是相当于要求由轴位置[x_left,x_right]描述的每个正方形的2个问题和[y_left,y_right]:

  1. 请告诉我的概率一个一维正态随机变量在区间[x_left,x_right]中?
  2. 请告诉我的另一概率,独立1D正态随机变量的存在在区间[y_left,y_right]?

因为这两个是独立的,为方形的全概率为:

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

这是一个确切的答案,所以计算时间不应该是一个问题!

编辑:我还要说,你可以选择在每个轴接近零更多蜱的网格,为了得到你想要的分辨率。以上每个广场的概率公式仍然成立。