如何使用 rpy2 在 Jupyter notebook 中应用 Henze-Zirkler 的多元正态性检验

Question

我有兴趣在 python 3x 中应用 Henze-Zirkler 的多元正态性检验，我想知道我是否可以在 Jupyter notebook 中的 python 中这样做。

我已经用我的数据拟合了一个 VAR 模型，然后我想测试这个拟合的 VAR 模型的残差是否服从正态分布。

如何使用 python 在 Jupyter notebook 中执行此操作？

Answer 1

R 中有一个包已经做了这个测试，它被称为 MVN

您要做的第一件事是将 MVN 导入 python，如

中所述

然后转到你的 jupyter notebook 并将 VAR(1) 模型拟合到你的数据中

# Fit VAR(1) Model

results = Model.fit(1)
results.summary()

将残差存储为 resi

resi=results.resid

然后

# Call function from R
import os
os.environ['R_USER'] = '...\Lib\site-packages\rpy2'
import rpy2.robjects as robjects
from rpy2.robjects import pandas2ri
pandas2ri.activate()

from rpy2.robjects.packages import importr

MVN = importr("MVN", lib_loc = "C:/.../R/win-library/3.3")

导入 MVN 后，您可以像这样简单地进行正态性测试

MVNresult =MVN.hzTest(resi, qqplot = 0)

如果你按

type(MVNresult)

你会发现它是一个

rpy2.robjects.methods.RS4

因此，在这种情况下，您会发现这个 link 在解释细节方面非常强大

之后

tuple(MVNresult.slotnames())

这将向您展示观察结果

('HZ', 'p.value', 'dname', 'dataframe')

那么你可能会得到这样的值

np.array(MVNresult.slots[tuple(MVNresult.slotnames())[i]])[0]

其中 i 代表 0, 1, 2, 3 为 'HZ', 'p-value',...

因此，如果 p 值即 i=1 小于 0.05，则残差 (resi) 在 5% 的置信水平下不是多元正态分布。

Answer 2

这是我后来发现的另一个答案。如果您不想将 R 的库导入 Python。可以将 R 的输出称为 python。也就是说，可以通过 python 激活 R 功能，如下所示：

import rpy2.robjects as robjects
from rpy2.robjects import r
from rpy2.robjects.numpy2ri import numpy2ri
from rpy2.robjects.packages import importr
import numpy as np

假设resi是python中的Dataframe say

# Create data
resi = pd.DataFrame(np.random.random((108, 2)), columns=['Number1','Number2'])

那么代码如下

#Converting the dataframe from python to R

# firt take the values of the dataframe to numpy
resi1=np.array(resi, dtype=float)

# Taking the variable from Python to R
r_resi = numpy2ri(resi1)

# Creating this variable in R (from python)
r.assign("resi", r_resi)

# Calling libraries in R 
r('library("MVN")')

# Calling a function in R (from python)
r("res <- hzTest(resi, qqplot = F)")

# Retrieving information from R to Python
r_result = r("res")

# Printing the output in python
print(r_result)

这将生成输出：

 Henze-Zirkler's Multivariate Normality Test 

--------------------------------------------- 

  data : resi 



  HZ      : 2.841424 

  p-value : 1.032563e-06 



  Result  : Data are not multivariate normal. 

---------------------------------------------

每 2021-08-25 更新 对 MVN 包和 ro rpy2 进行了一些 API 更改。以下适用于 MVN 5.9 版和 rpy2 3.4 版。

"""Interface file to access the R MVN package"""

import numpy as np
import rpy2.robjects.packages as rpackages
from rpy2.robjects import numpy2ri
from rpy2.robjects.packages import importr
from rpy2.robjects.vectors import StrVector


# Install packages, if they are not already installed
packages_to_install_if_needed = ("MVN",)
utils = rpackages.importr("utils")
utils.chooseCRANmirror(ind=1)  # select the first mirror in the list
names_to_install = [x for x in packages_to_install_if_needed if not rpackages.isinstalled(x)]
if len(names_to_install) > 0:
    utils.install_packages(StrVector(names_to_install))

# load the package
mvn = importr("MVN")

# Generate data
np_arr = np.random.multivariate_normal(np.ones(2), np.eye(2), size=100)

# activate automatic conversion from numpy to rpy2 interface objects
numpy2ri.activate()

# perform the work
res = mvn.mvn(np_arr)
print(res)

正在输出

$multivariateNormality
           Test        HZ   p value MVN
1 Henze-Zirkler 0.3885607 0.8343017 YES

$univariateNormality
              Test  Variable Statistic   p value Normality
1 Anderson-Darling  Column1     0.2443    0.7569    YES
2 Anderson-Darling  Column2     0.3935    0.3692    YES

$Descriptives
    n      Mean   Std.Dev    Median       Min      Max      25th     75th
1 100 0.9619135 1.0353688 1.0222279 -1.994833 3.679615 0.2696537 1.758255
2 100 0.7664778 0.9134449 0.8121996 -1.568635 2.648268 0.2068718 1.418113
        Skew    Kurtosis
1 -0.2123274 -0.16171832
2 -0.3718904 -0.05279222

Answer 3

有一个名为 Pingouin 的开源 Python 软件包，它提供 Henze-Zirkler 多元正态性检验并针对 R 的 MVM 进行了检验。

https://pingouin-stats.org/generated/pingouin.multivariate_normality.html

从文档中提取的示例：

import pingouin as pg
data = pg.read_dataset('multivariate')
X = data[['Fever', 'Pressure', 'Aches']]
pg.multivariate_normality(X, alpha=.05)
>>> HZResults(hz=0.5400861018514641, pval=0.7173686509624891, normal=True)

如何使用 rpy2 在 Jupyter notebook 中应用 Henze-Zirkler 的多元正态性检验

How to apply Henze-Zirkler's Multivariate Normality Test in Jupyter notebook with rpy2

python-3.x

statistics

jupyter-notebook

multivariate-testing