Я хочу создать гипотетические данные на основе корреляционной матрицы

Question

Я пытаюсь написать функцию, которая создает гипотетические данные. Функция получает корреляционную матрицу (R) и целое число (n) в качестве входных значений. Он создаст фрейм данных (или матрицу) с n строками, и двумерные корреляции между столбцами должны быть такими же, как в матрице корреляции (R).

Вот что у меня есть на данный момент:

Моя функция:

generate_data = function(R, n) {
r.eigen = eigen(R)
factors = t(t(r.eigen$vectors) * sqrt(r.eigen$values))
data = matrix(rnorm(n * ncol(R)), n)
data = data %*% t(factors)
return(data)}

И вот результаты, которые я получаю с разными n.

Пример матрицы корреляции:

R = matrix(c(1, .06, -.1, .1, .06, 1, -.51, .14, -.1, -.51, 1, .12, .1,  .14, .12, 1), ncol = 4)

> R
      [,1]  [,2]  [,3] [,4]

[1,]  1.00  0.06 -0.10 0.10

[2,]  0.06  1.00 -0.51 0.14

[3,] -0.10 -0.51  1.00 0.12

[4,]  0.10  0.14  0.12 1.00

А вот матрицы корреляции, которые я могу получить на основе n (количества строк).

>  round(cor(generate_data(R, 100)), 2)

      [,1]  [,2]  [,3] [,4]

[1,]  1.00 -0.23  0.09 0.12

[2,] -0.23  1.00 -0.44 0.23

[3,]  0.09 -0.44  1.00 0.09

[4,]  0.12  0.23  0.09 1.00

 >  round(cor(generate_data(R, 1000)), 2)

      [,1]  [,2]  [,3] [,4]

[1,]  1.00  0.05 -0.11 0.10

[2,]  0.05  1.00 -0.51 0.13

[3,] -0.11 -0.51  1.00 0.17

[4,]  0.10  0.13  0.17 1.00

 >  round(cor(generate_data(R, 10000)), 2)

      [,1]  [,2]  [,3] [,4]

[1,]  1.00  0.05 -0.09 0.10

[2,]  0.05  1.00 -0.50 0.13

[3,] -0.09 -0.50  1.00 0.14

[4,]  0.10  0.13  0.14 1.00

Моя функция, кажется, достаточно хорошо работает для больших n, но не работает для малых n. Есть ли способ создать функцию, которая работает и для меньших n?

Надеюсь, это достаточно ясно. Я ценю любую помощь.

Answer 1

Вы можете использовать функцию rmvnorm() из пакета mvtnorm. Однако для небольших n вариация выборки вряд ли стабилизируется, и вы получите маленькую вариацию выборки, которую вы также видели с помощью своей собственной функции. С этим ничего не поделаешь - это случайность.

library("mvtnorm")
R <- matrix(c(1, .06, -.1, .1, .06, 1, -.51, .14, -.1, -.51, 1, .12, .1,  .14, .12, 1), ncol = 4)
x <- rmvnorm(n=500, mean=c(0,0,0,0), sigma=R)

Это дает

cor(x)
           [,1]       [,2]        [,3]       [,4]
[1,]  1.0000000  0.1023989 -0.10946186 0.12230412
[2,]  0.1023989  1.0000000 -0.53853097 0.15985618
[3,] -0.1094619 -0.5385310  1.00000000 0.05587178
[4,]  0.1223041  0.1598562  0.05587178 1.00000000

, что кажется довольно близким к тому, что вы указали с помощью R. Если я увеличу количество образцов до 10000, я получу

x <- rmvnorm(n=10000, mean=c(0,0,0,0), sigma=R)
cor(x)
            [,1]        [,2]        [,3]      [,4]
[1,]  1.00000000  0.05969971 -0.08121426 0.1121826
[2,]  0.05969971  1.00000000 -0.51305601 0.1247779
[3,] -0.08121426 -0.51305601  1.00000000 0.1340828
[4,]  0.11218257  0.12477793  0.13408277 1.0000000