【原创】用python做Permutation Test置换检验

Posted 2023-05-11

参考技术A 一、概况

基本概念：Permutation test 置换检验是Fisher于20世纪30年代提出的一种基于大量计算（computationally intensive），利用样本数据的全（或随机）排列，进行统计推断的方法。

优势在于小样本检验：研究表明，当样本含量较大时, Permutation test得到的结果与经典的参数检验(t 检验、F 检验)近似。当样本含量较小时，Permutation test要优于参数检验，并且其检验效能也高于秩和检验。

原理：在具体使用上它和Bootstrap Methods类似，通过对样本进行顺序上的置换，重新计算统计检验量，构造经验分布，然后在此基础上求出P-value进行推断。

二、 实例

实验目的：验证加入某种生长素后拟南芥的侧根数量会明显增加。

实验设计：A组是加入某种生长素后，拟南芥的侧根数量；B是不加生长素时，拟南芥的侧根数量（均为假定值）。

  A组侧根数量（共12个数据）：24 43 58 67 61 44 67 49 59 52 62 50

  B组侧根数量（共16个数据）：42 43 65 26 33 41 19 54 42 20 17 60 37 42 55 28

检验方法：我们来用假设检验的方法来判断生长素是否起作用。

    我们的零假设H0为：加入的生长素不会促进拟南芥的根系发育。在这个检验中，若H0成立，那么A组数据的分布和B组数据的分布是一样的，A组数据和B组数据不存在显著差异，也就是服从同个分布。

    接下来构造检验统计量——A组侧根数目的均值同B组侧根数目的均值之差。

    statistic:= mean(Xa)-mean(Xb)

    对于观测值有 Sobs:=mean(Xa)-mean(Xb)=      (24+43+58+67+61+44+67+49+59+52+62+50)/12-  (42+43+65+26+33+41+19+54+42+20+17+60+37+42+55+28)/16=14   

    我们可以通过Sobs在置换分布（permutation distribution）中的位置来得到它的P-value。

    如果p<0.05，那么说明在原假设成立的情况下，出现这个Sobs值的概率是很低的(往极限讲的话在原假设成立的情况下是不会出现这个sob值的，那么既然现在这个值出现的，就可以拒绝原假设)，因此拒绝原假设，认为A、B两组数据存在显著差异，因此加入生长素会促进拟南芥的根系发育；

    如果p>0.05，那么说明在原假设成立的情况下，出现这个Sobs值的概率很大，原假设成立，认为A、B两组数据不存在显著差异，因此加入生长素不会促进拟南芥的根系发育。

检验过程：

Permutation test的具体步骤是：

1.将A、B两组数据合并到一个集合中，从中挑选出12个作为A组的数据（X'a），剩下的作为B组的数据（X'b）。

Group:=24 43 58 67 61 44 67 49 59 52 62 50 42 43 65 26 33 41 19 54 42 20 17 60 37 42 55 28

挑选出 X'a:=43 17 44 62 60 26 28 61 50 43 33 19

X'b:=55 41 42 65 59 24 54 52 42 49 37 67 67 20 42 58

2.计算并记录第一步中A组同B组的均值之差。Sper:=mean(X'a)-mean(X'b)= -7.875

3.对前两步重复999次（重复次数越多，得到的背景分布越”稳定“）

这样我们得到有999个置换排列求得的999个Sper结果，这999个Sper结果能代表拟南芥小样本实验的抽样总体情况。

import numpy as np

import matplotlib.pyplot as plt

import seaborn as sns

def exact_mc_perm_test(xs, ys, nmc):

    n, k = len(xs), 0

    diff = np.abs(np.mean(xs) - np.mean(ys))

    zs = np.concatenate([xs, ys])

    list=np.empty(nmc)

    for j in range(999):

        np.random.shuffle(zs)

        list[j]=np.abs(np.mean(zs[:n]) - np.mean(zs[n:]))

        k += diff < np.abs(np.mean(zs[:n]) - np.mean(zs[n:]))

    return list

xs = np.array([24,43,58,67,61,44,67,49,59,52,62,50])

ys = np.array([42,43,65,26,33,41,19,54,42,20,17,60,37,42,55,28])

list_a=exact_mc_perm_test(xs, ys, 999)

print(list_a)

sns.set_palette("hls") #设置所有图的颜色，使用hls色彩空间

sns.distplot(list_a,color="r",bins=30,kde=True) #kde=true，显示拟合曲线

plt.title('Permutation Test')

plt.xlabel('difference')

plt.ylabel('distribution')

plt.show()

如上图所示，我们的观测值 Sobs=14 在抽样总体右尾附近，说明在零假设条件下这个数值是很少出现的。在permutation得到的抽样总体中大于14的数值有9个，所以估计的P-value是9/999=0.01

最后还可以进一步精确P-value结果（做一个抽样总体校正），在抽样总体中加入一个远大于观测值 Sobs=14的样本，最终的P-value=(9+1)/(999+1)=0.01

结果表明我们的原假设不成立，加入生长素起到了促使拟南芥的根系发育的作用。

参考文献：

http://www.iikx.com/news/statistics/1824.html

https://www.plob.org/article/3176.html

尝试使用 next_permutation 在 C++ 中模拟 python 组合

【中文标题】尝试使用 next_permutation 在 C++ 中模拟 python 组合

【英文标题】：Trying to simulate python combinations in C++ with next_permutation

【发布时间】：2012-11-05 13:54:21

【问题描述】：

我需要将一个用 Python 编写的 sn-p 移植到 C++ 但是那个 sn-p 使用的是 python 中 itertools 的组合。

我真正有兴趣移植到 C++ 的那一行是：

for k in combinations(range(n-i),2*i):

Python 中的range(n-i) 将从0 to (n-i) - 1 生成一个列表

设 n = 16, i = 5

print range(n-i)

输出：

[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

python 组合将在该列表中生成所有可能的组合。

例如

print list(combinations(range(n-i),2*i))

输出：

[(0, 1, 2, 3, 4, 5, 6, 7, 8, 9), (0, 1, 2, 3, 4, 5, 6, 7, 8, 10), (0, 1, 2, 3, 4, 5, 6, 7, 9, 10), (0, 1, 2, 3, 4, 5, 6, 8, 9, 10), (0, 1, 2, 3, 4, 5, 7, 8, 9, 10), (0, 1, 2, 3, 4, 6, 7, 8, 9, 10), (0, 1, 2, 3, 5, 6, 7, 8, 9, 10), (0, 1, 2, 4, 5, 6, 7, 8, 9, 10), (0, 1, 3, 4, 5, 6, 7, 8, 9, 10), (0, 2, 3, 4, 5, 6, 7, 8, 9, 10), (1, 2, 3, 4, 5, 6, 7, 8, 9, 10)]

我想在 C++ 中使用 std::vector 和 next_permutation 生成类似的输出，但我仍然得到错误的结果。这是我目前的做法：

for(int j = 0; j < n-i; j++) 
        temp_vector.push_back(j); 

那个 sn-p 相当于 Python 中的range(n-i)。

但是下面的sn-p：

do 
     myvector.push_back(temp_vector);
 while(next_permutation(temp_vector.begin(),temp_vector.begin()+2*i));
cout<<myvector.size()<<endl;

不等同于 Python 中的combinations(range(n-i),2*i))，我尝试了很多变体，但仍然无法得出我期望的结果。

例如：

设 n = 16 i = 5

Python

>>> print len(list(combinations(range(n-i),2*i)))

11

C++

#include <vector>
#include <iostream>

using namespace std;
int main() 
    vector<int> temp_vector;
    vector< vector<int> > myvector;
    int n = 16, i = 5;
    for(int j = 0; j < n - i; j++) 
            temp_vector.push_back(j);
    
    do 
            myvector.push_back(temp_vector);
     while(next_permutation(temp_vector.begin(), temp_vector.begin()+2*i));
    cout<<myvector.size()<<endl;
    return 0;

g++ combinations.cpp

./a.out

3628800

任何指导将不胜感激！非常感谢！

【参考方案1】：

组合和排列不是一回事。

组合是另一个集合中项目子集的无序列表。排列是列表中项目的唯一顺序。

您正在从 11 项内容的列表中生成 10 项内容的所有组合，因此您将获得 11 个结果，每个结果都缺少原始 11 项中的一项。

生成每个排列将生成原始 11 项的每个唯一顺序。由于本例中的项目都是唯一的，这意味着结果将是 11！列出每个包含所有 11 个项目的位置。但是，您仅从前 10 个项目生成排列，因此您将获得 10 个！列表，其中不包含第 11 项。

您需要找到一种算法来生成组合而不是排列。

---

没有用于组合的内置算法。 std::next_permutation 可用作生成组合的算法的一部分：参见Generating combinations in c++。

Here's 组合算法的旧提案草案，包括代码。

【讨论】：

我查看代码已经很久了，以至于我没有意识到 C++ 算法正在生成排列，而不是组合......我现在觉得很愚蠢！谢谢你让我大开眼界！

您是否碰巧知道在 C++ 中生成组合的内置算法？这对我有很大帮助:)