2015-07-04 129 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)) 

這是一個確切的答案,所以計算時間不應該是一個問題!

編輯:我還要說,你可以選擇在每個軸接近零更多蜱的網格,爲了得到你想要的分辨率。以上每個廣場的概率公式仍然成立。