Как сгенерировать многомерное распределение в J, учитывая вектор среднего значения и ковариационную матрицу?
Например, в Python np.random.multivariate_normal([0,0],[[1,.75],[.75,1]],1000) генерирует многомерное распределение с [0,0] в качестве вектора среднего значения и [[1,.75],[.75,1]] в качестве матрицы дисперсии-ковариации? Спасибо
Википедия описывает стандартный метод создания многомерных нормальных распределений:
Mu=: 0 0 NB. vector of means
]Sigma=: 1 0.75 ,: 0.75 1 NB. covariance matrix
1 0.75
0.75 1
Мы можем получить разложение Холецкого ковариационной матрицы, используя код из скрипта matfacto.ijs в аддоне math/misc (или использовать аддон LAPACK)
load 'math/misc/matfacto'
A=: choleski Sigma NB. Cholesky decomp
Создайте 2 независимые одномерные нормальные переменные, используя дополнение stats/distribs.
load 'stats/distribs'
z=: rnorm 2 1000 NB. 2 standard normal variables sampled 1000 times
Теперь сгенерируйте желаемые многомерные распределения:
X=: Mu + A mp z
Теперь проверьте, что дистрибутивы соответствуют указанным:
load 'stats/base'
mean"1 X
0.0264368 0.00887907 NB. mean close to 0 (Mu)
stddev"1 X
0.987214 0.991614 NB. stddev close to 1 (sqrt of diagonal of Sigma)
corr/ X
0.746917 NB. correlation close to 0.75 (off-diagonal of Sigma)
Мы можем закодировать это как один глагол:
multivar_norm=: dyad define
'Mu Sigma'=. x
A=. choleski Sigma
z=. rnorm y ,~ #Sigma
Mu + A mp z
)
X=: (Mu;Sigma) multivar_norm 1000
((mean , stddev)"1 ; corr/) X
┌──────────────────┬────────┐
│0.0199138 1.01788│0.749184│
│ 0.035176 0.987191│ │
└──────────────────┴────────┘
