2011-04-15 37 views
5

我目前使用python和RPY裏面使用R.功能基於R兌現爲LatinHypercube /蒙特卡洛試驗相關

如何使用[R庫生成兌現2個變量之間的相關性蒙特卡羅樣本.. 例如 如果變量A和B具有85%(0.85)的相關性,我需要生成所有的蒙特卡洛樣品履行之間的& B.

希望如果任何人可以共享想法/片段的相關性

謝謝

+0

我很高興你提出這個問題,因爲我必須爲一個項目做一些敏感性分析,並且一直在研究這個理論,但還沒有搜索出R選項來實現我正在學習的方法。 – 2011-04-15 16:24:25

回答

5

這是一個FAQ。下面是一個使用推薦的包一個答案:

R> library(MASS) 
R> example(mvrnorm) 

mvrnrmR> Sigma <- matrix(c(10,3,3,2),2,2) 

mvrnrmR> Sigma 
    [,1] [,2] 
[1,] 10 3 
[2,] 3 2 

mvrnrmR> var(mvrnorm(n=1000, rep(0, 2), Sigma)) 
     [,1] [,2] 
[1,] 8.82287 2.63987 
[2,] 2.63987 1.93637 

mvrnrmR> var(mvrnorm(n=1000, rep(0, 2), Sigma, empirical = TRUE)) 
    [,1] [,2] 
[1,] 10 3 
[2,] 3 2 
R> 

相關和協方差之間的切換是直接的(提示:標準偏差的矢量的外積)。

7

Iman and Conover的排名相關方法似乎是一個廣泛使用和一般的方法來產生相關蒙特卡洛樣本進行基於計算機的實驗,敏感性分析等。不幸的是,我只是剛剛遇到這個,並沒有訪問PDF所以不知道作者實際上如何實施他們的方法,但你可以按照這一點。

他們的方法更通用,因爲每個變量都可以來自不同的分佈,而不像@Dirk答案的多元正態分佈。

更新:我在包mc2d中發現了上述方法的R實現,特別是你想要的cornode()函數。

下面是從?cornode

> require(mc2d) 
> x1 <- rnorm(1000) 
> x2 <- rnorm(1000) 
> x3 <- rnorm(1000) 
> mat <- cbind(x1, x2, x3) 
> ## Target 
> (corr <- matrix(c(1, 0.5, 0.2, 0.5, 1, 0.2, 0.2, 0.2, 1), ncol=3)) 
    [,1] [,2] [,3] 
[1,] 1.0 0.5 0.2 
[2,] 0.5 1.0 0.2 
[3,] 0.2 0.2 1.0 
> ## Before 
> cor(mat, method="spearman") 
      x1   x2   x3 
x1 1.00000000 0.01218894 -0.02203357 
x2 0.01218894 1.00000000 0.02298695 
x3 -0.02203357 0.02298695 1.00000000 
> matc <- cornode(mat, target=corr, result=TRUE) 
Spearman Rank Correlation Post Function 
      x1  x2  x3 
x1 1.0000000 0.4515535 0.1739153 
x2 0.4515535 1.0000000 0.1646381 
x3 0.1739153 0.1646381 1.0000000 

採取matc的等級相關現在非常接近的corr目標的相關性的示例。

這樣做的想法是,您從每個變量的分佈中分別抽取樣本,然後使用Iman & Connover方法儘可能使樣本(儘可能接近)與目標相關性相近。

+0

好的答案,加文,因爲我們從來沒有被告知這是否只是多元正態(或t或...),或者,正如你猜測的,有點兒聰明的東西,這是一個不明確的問題。 – 2011-04-15 18:01:53

+0

@Dirk謝謝 - 我最初的想法是迴應你對`mvrnorm()`所做的反應,但是你打敗了我,然後我想知道這個我一直在閱讀的等級相關方法。 – 2011-04-15 18:05:43

2

此問題未被標記爲python,但根據您的評論,它看起來像你也可能正在尋找一個Python解決方案。最基本的Python實現伊曼Convover的,我可以編造看起來像在Python以下(實際上numpy的):

def makeCorrelated(y, corMatrix): 
    c = multivariate_normal(zeros(size(y, 0)) , corMatrix, size(y, 1)) 
    key = argsort(argsort(c, axis=0), axis=0).T 
    out = map(take, map(sort, y), key) 
    out = array(out) 
    return out 

其中y是樣品的從邊緣分佈的陣列和corMatrix是正半定,對稱相關矩陣。鑑於這個函數對c矩陣使用了multivariate_normal(),你可以告訴它使用隱含的高斯Copula。要使用不同的copula結構,您需要爲c矩陣使用不同的驅動程序。