random - 如何在 J 中生成多元正态分布

Question

给定均值向量和协方差矩阵，如何在 J 中生成多元分布？

例如，在 Python 中， np.random.multivariate_normal([0,0],[[1,.75],[.75,1]],1000) 生成一个多元分布，其中 [0,0] 作为平均值向量，[[1,.75],[.75,1]] 作为方差-协方差矩阵？ 谢谢

Answer 1

维基百科描述了 standard method for creating multivariate normal distributions :

   Mu=: 0 0                   NB. vector of means
   ]Sigma=: 1 0.75 ,: 0.75 1  NB. covariance matrix
   1 0.75
0.75    1

我们可以使用 math/misc 插件中的 matfacto.ijs 脚本中的代码获得协方差矩阵的 Cholesky 分解(或使用 LAPACK addon )

   load 'math/misc/matfacto'
   A=: choleski Sigma         NB. Cholesky decomp

使用 stats/distribs 创建 2 个独立的单变量正态变量插件。

   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│        │
└──────────────────┴────────┘