Генерация пар случайных чисел, равномерно распределенных и коррелированных

14

Я хотел бы генерировать пары случайных чисел с определенной корреляцией. Однако обычный подход использования линейной комбинации двух нормальных переменных здесь недопустим, поскольку линейная комбинация равномерных переменных больше не является равномерно распределенной переменной. Мне нужно, чтобы две переменные были одинаковыми.

Любая идея о том, как генерировать пары однородных переменных с заданной корреляцией?

Onturenio
источник
6
Близко связаны: stats.stackexchange.com/questions/30526 . Вы также хотите проверить тег связки - просто нажмите на ссылку здесь. Быстрая и грязная техника состоит в том, чтобы X был равномерным [0,1] и Y=X когда Xα и Y=1+αX противном случае. Соотношение ρ=2(α1)3+1 , откуда α=1((1ρ)/2)1/3 делает трюк. Но связки дадут вам больше контроля ...
whuber
Спасибо за комментарий, но да, я думаю, что этот метод действительно "грязный"
Onturenio
1
Я надеялся, что, увидев этот подход, вы поймете, что можете (и должны) предоставить дополнительные критерии, касающиеся свойств ваших пар случайных чисел. Если это «грязно», то что именно не так с решением? Скажите нам, чтобы мы могли предоставить более подходящие ответы для вашей ситуации.
whuber
На этот вопрос был дан ответ случайно в ответ на тесно связанный вопрос: как генерировать пары RV с линейной регрессионной зависимостью. Поскольку наклон линейной регрессии легко вычисляется с коэффициентом корреляции, и все возможные наклоны могут быть получены, это дает способ произвести именно то, что вы хотите. См. Stats.stackexchange.com/questions/257779/… .
whuber
1
Также см. Stats.stackexchange.com/questions/31771 , в котором дано обобщение на три случайные формы.
whuber

Ответы:

16

Я не знаю универсального метода для генерации коррелированных случайных величин с любым заданным предельным распределением. Итак, я предложу специальный метод для генерации пар равномерно распределенных случайных величин с заданной (Pearson) корреляцией. Без ограничения общности я предполагаю, что желаемое предельное распределение является стандартным равномерным (т. Е. Поддержка [0,1] ).

Предлагаемый подход основан на следующем:
a) Для стандартных равномерных случайных величин и U 2 с соответствующими функциями распределения F 1 и F 2 имеем F i ( U i ) = U i , для i = 1 , 2 . Таким образом, по определению число Спирмена равно ρ S ( U 1 , U 2 ) = c o r r ( FU1U2F1F2Fi(Ui)=Uii=1,2 Таким образом, коэффициенты корреляции Спирмена и Пирсона равны (примерные версии могут отличаться).

ρS(U1,U2)=corr(F1(U1),F2(U2))=corr(U1,U2).

б) Если являются случайными величинами с непрерывными полями и гауссовой копулой с коэффициентом корреляции (Пирсона) ρ , то число Спирмена равно ρ S ( X 1 , X 2 ) = 6X1,X2ρ Это позволяет легко генерировать случайные величины, которые имеют желаемое значение ро Спирмена.

ρS(X1,X2)=6πarcsin(ρ2).

Подход заключается в том, чтобы генерировать данные из гауссовой связки с подходящим коэффициентом корреляции , так что относительное число Спирмена соответствует желаемой корреляции для однородных случайных величин.ρ

Алгоритм моделирования
Пусть обозначает желаемый уровень корреляции, а n - количество генерируемых пар. Алгоритм:rn

  1. Вычислить .ρ=2sin(rπ/6)
  2. Генерация пары случайных величин из гауссовой связки (например, при таком подходе )
  3. Повторите шаг 2 раз.n

Пример
Следующий код является примером реализации этого алгоритма с использованием R с целевой корреляцией и n = 500 пар.r=0.6n=500

## Initialization and parameters 
set.seed(123)
r <- 0.6                            # Target (Spearman) correlation
n <- 500                            # Number of samples

## Functions
gen.gauss.cop <- function(r, n){
    rho <- 2 * sin(r * pi/6)        # Pearson correlation
    P <- toeplitz(c(1, rho))        # Correlation matrix
    d <- nrow(P)                    # Dimension
    ## Generate sample
    U <- pnorm(matrix(rnorm(n*d), ncol = d) %*% chol(P))
    return(U)
}

## Data generation and visualization
U <- gen.gauss.cop(r = r, n = n)
pairs(U, diag.panel = function(x){
          h <- hist(x, plot = FALSE)
          rect(head(h$breaks, -1), 0, tail(h$breaks, -1), h$counts/max(h$counts))})

На рисунке ниже, диагональные графики показывают гистограммы переменных и U 2 , а недиагональные графики показывают графики рассеяния U 1 и U 2 . U1U2U1U2введите описание изображения здесь

По построению случайные величины имеют одинаковые поля и коэффициент корреляции (близкий к) . Но из-за эффекта выборки коэффициент корреляции смоделированных данных не точно равен r .rr

cor(U)[1, 2]
# [1] 0.5337697

Обратите внимание, что gen.gauss.copфункция должна работать с более чем двумя переменными, просто указав большую корреляционную матрицу.


r=0.5,0.1,0.6n

## Simulation
set.seed(921)
r <- 0.6                                                # Target correlation
n <- c(10, 50, 100, 500, 1000, 5000); names(n) <- n     # Number of samples
S <- 1000                                               # Number of simulations

res <- sapply(n,
              function(n, r, S){
                   replicate(S, cor(gen.gauss.cop(r, n))[1, 2])
               }, 
               r = r, S = S)
boxplot(res, xlab = "Sample size", ylab = "Correlation")
abline(h = r, col = "red")

введите описание изображения здесь введите описание изображения здесь введите описание изображения здесь

QuantIbex
источник
3
Общий метод генерации коррелированных многомерных распределений с заданными маргинальными распределениями называется связкой .
whuber
@whuber, использование связок позволяет определить структуру зависимости между случайными переменными. Проблема в том, что на (человека) корреляцию влияют как структура зависимости, так и поля. Таким образом, каждый выбор полей потребует соответствующего выбора параметров связки, не говоря уже о том, что некоторые уровни корреляции просто не могут быть достигнуты для данных полей (например, см. Здесь ). Если вам известен метод, который позволяет «контролировать» уровень корреляции для любого выбора полей, я хотел бы узнать об этом.
QuantIbex
Спасибо @QuantIbex. Но я не понимаю, почему «а) подразумевает, что коэффициенты корреляции Ро и (Пирсона) Спирмена для случайных величин со стандартными однородными полями приблизительно равны в большой выборке»
Онтуренио,
2
[1,1]
1
@Quantibex Я позволил себе добавить предложение, которое указывает, что ваша gen.gauss.copфункция будет работать для более чем двух переменных с (тривиальной) настройкой. Если вам не нравится дополнение или вы хотите поставить его по-другому, пожалуйста, отмените или измените его по мере необходимости.
Glen_b
0

u1U(0,1)u1w1U(0,1)I=1u1w2U(0,1)I=0u1U(0,1)u2

E(u1u2)=E[Iw1+(1I)w2][Iw1+(1I)w3]

I(I1)=0I2=I(1I)2=(1I)I01Iw

E(u1u2)=E(I)E(w12)+E(1I)E(w2)E(w3) =pE(w12)+(1p)/4

V(w1)=1/12E(w12)=1/3E(u1u2)=p/12+1/4cov(u1u2)=p/12V(u1)=V(u2)=1/12cor(u1,u2)=p

Нил Оден
источник
0

(u1,u2)=Iw1+(1I)(w2,w3), where w1,w2, and w3 are independent U(0,1) and I is Bernoulli(p). u1 and u2 will then have U(0,1) distributions with correlation p. This extends immediately to k-tuples of uniforms with compound symmetric variance matrix.

If you want pairs with negative correlation, use (u1,u2)=I(w1,1w1)+(1I)(w2,w3), and the correlation will be p.

Neal Oden
источник
Can you add a short proof of why this works?
The Laconic
if your want to be computationally efficient, u1=w1 also produces the same correlation (both positive and negative cases)
Anvit