Как вызвать корреляции между двумя обратными кумулятивными распределениями вероятностей в [r]? - PullRequest
2 голосов
/ 18 марта 2020

Я хотел бы создать коррелированное обратное совокупное распределение. В настоящее время, например, у меня есть два обратных распределения, показанных ниже, но я хотел бы, например, вызвать корреляцию, скажем, -0,5. Есть ли способ, которым я могу достичь этого?


library(lognorm)
library(dplyr)

Var_a <- tbl_df(qlnorm(runif(1000), meanlog = 0.0326, sdlog = 0.0288))
var_b <- tbl_df(qlnorm(runif(1000), meanlog = 0.0452, sdlog = 0.0364))

cor(Var_a, var_b)

Ответы [ 2 ]

1 голос
/ 18 марта 2020

Если у вас есть 15 переменных с матрицей корреляции CC, вы можете использовать гауссову связку, чтобы получить коррелированные равномерные переменные, используя разложение Холецкого CC, а затем инвертировать те, у которых указаны указанные маргиналы, как вы делали выше. ( См. Здесь , например).

nv <- NROW(CC)
num_samples <- 1000
A <- matrix(rnorm(num_samples * nv), ncol = nv)
U <- pnorm(A %*% chol(CC))

Если ваши 15 переменных имеют свои средние значения и стандартные отклонения, сохраненные в векторах means и stdevs, вы можете сделать:

rv <- sapply(1:nv, function(i) qlnorm(U[,i], meanlog = means[i], sdlog = stdevs[i]))

rv - это ваши смоделированные переменные с близкой к желаемой структурой корреляции, которую вы можете проверить с помощью cor(rv).

1 голос
/ 18 марта 2020

Будет ли работать следующее для вас?

set.seed(100)
x1 <- rnorm(1000)
y1 <- rnorm(1000) - .6 * x1

x2 = pnorm(x1)
y2 = pnorm(y1)

cor(cbind(x2, y2))
#            x2         y2
# x2  1.0000000 -0.4995593
# y2 -0.4995593  1.0000000

Var_a <- tbl_df(qlnorm(x2, meanlog = 0.0326, sdlog = 0.0288))
var_b <- tbl_df(qlnorm(y2, meanlog = 0.0452, sdlog = 0.0364))

cor(Var_a, var_b)
#            value
# value -0.5239145

update: все еще не понимаете, что вы делаете, но если вы просто хотите применить то, что я сделал, к 15 переменным, возможно, сделайте что-то подобное?

library(MASS)
sigma <- matrix(.5, nrow = 15, ncol = 15) + diag(15)*.5  #your correlation matrix
sigma
vars <- mvrnorm(1000, mu = rep(0, 15), Sigma = sigma)
vars
cor(vars)
vars2 <- pnorm(vars)
cor(vars2)
#use each of these as variable in qlnorm

vars2 <- data.frame(vars2)
names(vars2)
vars2

vars2[paste("log_", 1:15)] <- lapply(vars2[, 1:15], function(x) {qlnorm(x, meanlog = 0.0326, sdlog = 0.0288)})
names(vars2)
vars2 <- vars2[, -c(1:15)]
cor(vars2)
...