数据准备
许多实际情况中统计假设（假定观测数据抽样自正态分布或者其他性质较好的理论分布）并不一定满足，比如数据抽样于未知或混合分布、样本量过小、存在离群点、基于理论分布设计合适的统计检验过于复杂且数学上难以处理等情况，这时基于随机化和重抽样的统计方法就可派上用场。
置换检验的定义
置换检验(Permutation test)，也称随机化检验或重随机化检验，是Fisher于20世纪30年代提出的一种基于大量计算（computationally intensive），利用样本数据的全（或随机）排列，进行统计推断的方法，因其对总体分布自由，应用较为广泛，特别适用于总体分布未知的小样本资料，以及某些难以用常规方法分析资料的假设检验问题。
置换检验的原理
1、提出原假设，比如XX处理后结果没有变化
2、计算统计量，如两组的均值之差，记作t0
3、将所有样本放在一起，然后随机排序进行分组，再计算其统计量t1
4、重复第3步骤，直至所有排序可能性都齐全（比如有A组有n样本，B组有m样本，则总重复次数相当于从n+m中随机抽取n个的次数），得到一系列的统计量（t1-tn）
5、最后将这些统计量按照从小到大排序，构成抽样分布，再看t0是否落在分布的置信区间内（如95%置信区间），这时候可计算一个P值（如果抽样总体1000次统计量中大于t0的有10个，则估计的P值为10/1000=0.01），落在置信区间外则拒绝原假设
6、如果第3步骤是将所有可能性都计算了的话，则是精确检验；如果只取了计算了部分组合，则是近似结果，这时一般用蒙特卡罗模拟（Monte Carlo simulation）的方法进行置换检验
7、置换检验和参数检验都计算了统计量，但是前者是跟置换观测数据后获得的经验分布进行比较，后者则是跟理论分布进行比较。
请牢记：置换检验都是使用伪随机数来从所有可能的排列组合中进行抽样（当做近似检验时）。因此，每次检验的结果都有所不同。
coin包提供了一个进行置换检验的一般性框架。通过该包，你可以回答如下问题。
响应值与组的分配独立吗？
两个数值变量独立吗？
两个类别型变量独立吗？
表12-2列出来的每个函数都是如下形式：
function_name(formula, data, distribution=)
其中：
formula描述的是要检验变量间的关系。示例可参见表12-2；
data是一个数据框；
distribution指定经验分布在零假设条件下的形式，可能值有exact，asymptotic和
approximate。
若distribution = "exact"，那么在零假设条件下，分布的计算是精确的（即依据所有可能的排列组合）。当然，也可以根据它的渐进分布（distribution = "asymptotic"）或蒙特卡洛重抽样（distribution = "approxiamate(B = #)"）来做近似计算，其中#指所需重复的次数。
distribution = "exact"当前仅可用于两样本问题。
传统t检验表明存在显著性差异（p <0.05），而精确检验却表明差异并不显著（p >0.072）。
第7章我用自己的数据进行了t检验，对比一下传统t检验和置换检验，结果如下：
两种检验方式下结果都是显著的
Wilcoxon-Mann-Whitney U检验
coin包规定所有的类别型变量都必须以因子形式编码。
wilcox.test()默认计算的也是精确分布。
K样本检验的置换检验
通过chisq_test()或cmh_test()函数，我们可用置换检验判断两类别型变量的独立性。 当数据可根据第三个类别型变量进行分层时，需要使用后一个函数。若变量都是有序型，可使用lbl_test()函数来检验是否存在线性趋势。
卡方独立性检验
卡方独立性检验的置换检验
你可能会有疑问，为什么需要把变量Improved从一个有序因子变成一个分类因子？（好问题！）这是因为，如果你用有序因子，coin()将会生成一个线性与线性趋势检验，而不是卡方检验。
结果解读：两种检验下p值都是小于0.05，说明Treatment和Improved之间相互不独立
自己数据的演示
结果解读：p值均为1，表明nitrogen和variety相互独立。
spearman_test()函数提供了两数值变量的独立性置换检验。
当处于不同组的观测已经被分配得当，或者使用了重复测量时，样本相关检验便可派上用场。
对于两配对组的置换检验，可使用wilcoxsign_test()函数；多于两组时，使用friedman_test()函数。
自己数据演示
lmPerm包可做线性模型的置换检验。比如lmp()和aovp()函数即lm()和aov()函数的修改版，能够进行置换检验，而非正态理论检验。
lmp()和aovp()函数的参数与lm()和aov()函数类似，只额外添加了perm =参数。
perm =选项的可选值有"Exact"、"Prob"或"SPR"。Exact根据所有可能的排列组合生成精确检验。Prob从所有可能的排列中不断抽样，直至估计的标准差在估计的p值0.1之下，判停准则由可选的Ca参数控制。SPR使用贯序概率比检验来判断何时停止抽样。注意，若观测数大于10，perm = "Exact"将自动默认转为perm = "Prob"，因为精确检验只适用于小样本问题。
简单线性回归的置换检验
R语言实战的例子：
多项式回归的置换检验
R语言实战的例子：
自己数据集的例子：
R语言实战的例子：
自己数据集的例子：
当两种方法所得结果不一致时，你需要更加谨慎地审视数据，这很可能是因为违反了正态性假设或者存在离群点。
R语言实战的例子：
自己数据集的例子：
R语言实战的例子：
自己数据集的例子：
R语言实战的例子：
自己数据集的例子：
值得注意的是，当将aovp()应用到方差分析设计中时，它默认使用唯一平方和法（SAS也称为类型III平方和）。每种效应都会依据其他效应做相应调整。R中默认的参数化方差分析设计使用的是序贯平方和（SAS是类型I平方和）。每种效应依据模型中先出现的效应做相应调整。对于平衡设计，两种方法结果相同，但是对于每个单元格观测数不同的不平衡设计，两种方法结果则不同。不平衡性越大，结果分歧越大。若在aovp()函数中设定seqs = TRUE，可以生成你想要的序贯平方和。
你可能已经注意到，基于正态理论的检验与上面置换检验的结果非常接近。在这些问题中数据表现非常好，两种方法结果的一致性也验证了正态理论方法适用于上述示例。
当然，置换检验真正发挥功用的地方是处理非正态数据（如分布偏倚很大）、存在离群点、样本很小或无法做参数检验等情况。不过，如果初始样本对感兴趣的总体情况代表性很差，即使是置换检验也无法提高推断效果。
置换检验主要用于生成检验零假设的p值，它有助于回答“效应是否存在”这样的问题。不过，置换方法对于获取置信区间和估计测量精度是比较困难的。幸运的是，这正是自助法大显神通的地方。
所谓自助法，即从初始样本重复随机替换抽样，生成一个或一系列待检验统计量的经验分布。 无需假设一个特定的理论分布，便可生成统计量的置信区间，并能检验统计假设。
倘若你假设均值的样本分布不是正态分布，该怎么办呢？可使用自助法。
(1)从样本中随机选择10个观测，抽样后再放回。有些观测可能会被选择多次，有些可能一直都不会被选中。
(2)计算并记录样本均值。
(3)重复1和2一千次。
(4)将1000个样本均值从小到大排序。
(5)找出样本均值2.5%和97.5%的分位点。此时即初始位置和最末位置的第25个数，它们就限定了95%的置信区间。
boot包扩展了自助法和重抽样的相关用途。你可以对一个统计量(如中位数)或一个统计量向量(如一列回归系数)使用自助法。
一般来说，自助法有三个主要步骤。
(1)写一个能返回待研究统计量值的函数。如果只有单个统计量（如中位数），函数应该返回一个数值；如果有一列统计量（如一列回归系数），函数应该返回一个向量。
(2)为生成R中自助法所需的有效统计量重复数，使用boot()函数对上面所写的函数进行处理。
(3)使用boot.ci()函数获取第(2)步生成的统计量的置信区间。
主要的自助法函数是boot()，它的格式为：
bootobject <- boot(data=, statistic=, R=, ...)
参数见下表：
boot()函数调用统计量函数R次，每次都从整数1:nrow(data)中生成一列有放回的随机指标，这些指标被统计量函数用来选择样本。统计量将根据所选样本进行计算，结果存储在bootobject中。
你可以用bootobject t0和bootobject t来获取这些元素。
一旦生成了自助样本，可通过print()和plot()来检查结果。如果结果看起来还算合理， 使用boot.ci()函数获取统计量的置信区间。格式如下:
boot.ci(bootobject, conf=, type= )
type参数设定了获取置信区间的方法。perc方法（分位数）展示的是样本均值，bca将根据偏差对区间做简单调整。
回归的R平方值
1000次自助抽样
输出结果
结果可视化
95%的置信区间获取
回归系数向量函数
自助抽样1000次
获得车重和发动机排量95%的置信区间
置换检验和自助法并不是万能的，它们无法将烂数据转化为好数据。当初始样本对于总体情况的代表性不佳，或者样本量过小而无法准确地反映总体情况，这些方法也是爱莫能助。
参考资料：
原理参考 文章 ，主要思想我认为是求出所有分布的可能（中间的一般为零假设），出现这种分布的概率。
distribution= 参数可为exact（精确模式，即依据所有可能的排列组合，仅适用于两样本问题）、approxiamate(nresample=#)（蒙特卡洛抽样，#指需要重复的次数）、asymptotic（渐进分布抽样）
lmPerm包更擅长方差分析。示例实验设计仍为5组接受不同治疗方法的多组结果比较。
实验示例仍为关节炎的治疗（两种）与效果（无、部分、显著）间的关系
实验示例为研究文盲率与谋杀率是否相关
主要为 lmp() 、 aovp() 两个函数分别对应参数法的 lm() 线性回归、 aov() 方差分析。主要格式上的区别是添加了 perm= 参数。可以为Exact（精确模式）、Prob（不断从可能的序列中抽样，直至估计的标准差在估计的p值0.1之下）、SPR（使用贯序概率比检验来判断何时停止抽样）。值得注意的是当样本观测大于10，perm="Exact"自动默认转为"Prob"，因此精确检验只适用于小样本问题。
（1）简单线性回归
实验示例仍为以身高预测体重的设计
（2）多项式回归
高精度拟合身高体重回归关系
（3）多元回归
探究谋杀率与多因素的回归关系
（1）单因素方差分析
（2）单因素协方差分析
实验示例仍为药物对刚出生小鼠体重影响，协变量为怀孕时间
（3）双因素方差分析（交互效应）
实验示例：两种药物分别在不同剂量下对小鼠牙齿长度的影响。
核心思想是有放回的抽样多次（1000次）
（1）写一个能返回带研究统计量的函数；
（2）确定重复数，使用 boot() 函数处理；（一般重复1000次即可；此外有人认为初始样本大小为20-30即可得到足够好的结果）；
（3） boot.ci() 函数计算统计量置信区间。
实验示例：使用mtcar数据框，采用多元回归，根据车重和发动机排量来预测汽车的每加仑行驶的英里数。想获得95%的R平方值（预测变量对响应变量可解释的方差比）的置信区间
（1）首先写函数
（2）然后使用boot()函数
（3）最后boot.ci()函数求置信区间
实验示例：使用mtcar数据框，采用多元回归，根据车总和发动机排量来预测汽车的每加仑行驶的英里数。想获取一个统计量向量--三个回归系数（截距项、车总、发动机排量）95%的置信区间。
--title: R语言中dnorm, pnorm, qnorm与rnorm以及随机数
date: 2018-09-07 12:02:00
type: "tags"
tags:
在R语言中，与正态分布（或者说其它分布）有关的函数有四个，分别为dnorm,pnorm,qnorm和rnorm，其中，dnorm表示密度函数，pnorm表示分布函数，qnorm表示分位数函数，rnorm表示生成随机数的函数。在R中与之类似的函数还有很多，具体的可以通过 help(Distributions) 命令去查看，对于分位数或百分位数的一些介绍可以看这篇笔记 《分位数及其应用》 ，关于正态分布的知识可以看这篇笔记 《正态分布笔记》 。
现在这篇笔记就介绍一下这些函数的区别。
R提供了多种随机数生成器(random number generators, RNG)，默认采用的是Mersenne twister方法产生的随机数，该方法是由Makoto Matsumoto和Takuji Nishimura于1997年提出来的，其循环周期是 。R里面还提供了了Wichmann-Hill、Marsaglia-Multicarry、Super-Duper、Knuth-TAOCP-2002、Knuth-TAOCP和L'Ecuyer-CMRG等几种随机数生成方法，可以通过 RNGkind() 函数进行更改，例如，如果要改为WIchmann-Hill方法，就使用如下语句：
在R中使用随机数函数，例如 rnorm() 函数来生成的随机数是不一样的，有时我们在做模拟时，为了比较不同的方法，就需要生成的随机数都一样，即重复生成相同的随机数，此时就可以使用 set.seed() 来设置随机数种子，其参数为整数，如下所示：
dnorm 中的 d 表示 density ， norm 表示正态贫，这个函数是正态分布的 概率密度(probability density)函数 。
正态分布的公式如下所示：
给定x，μ和σ后， dnorm() 这个函数返回的就是会返回上面的这个公式的值，这个值就是Z-score，如果是标准正态分布，那么上述的公式就变成了这个样子，如下所示：
现在看一个案例，如下所示：
dnorm(0,mean=0,sd=1) 由于是标准正态分布函数的概率密度，这个命令其实可以直接写为 dnorm(0) 即可，如下所示：
再看一个非标准正态分布的案例，如下所示：
虽然在 dnorm() 中，x是一个概率密度函数(PDF,Probability Density Function)的独立变量，但它也能看作是一组经过Z转换后的一组变量，现在我们看一下使用 dnorm 来绘制一个正态分布的概率密度函数曲线，如下所示：
现在使用 dnorm() 函数计算一下Z_scores的概率密度，如下所示：
现在绘图，如下所示：
从上面的结果可以看出，在每个Z-score处， dnorm 可以绘制出这个Z-score对应的正态分布的pdf的高度。
pnorm 函数中的 p 表示Probability，它的功能是，在正态分布的PDF曲线上，返回从负无穷到 q 的积分，其中这个 q 指的是一个Z-score。现在我们大概就可以猜测出 pnorm(0) 的值是0.5，因为在标准正态分布曲线上，当Z-score等于0时，这个点正好在标准正态分布曲线的正中间，那么从负无穷到0之间的曲线面积就是整个标准正态分布曲线下面积的一半，如下所示：
pnorm 函数还能使用 lower.tail 参数，如果 lower.tail 设置为 FALSE ，那么 pnorm() 函数返回的积分就是从 q 到正无穷区间的PDF下的曲线面积，因此我们就知道了， pnorm(q) 与 1-pnorm(q,lower.tail=FALSE) 的结果是一样的，如下所示：
在计算机出现之前的时代里，统计学家们使用正态分布进行统计时，通常是要查正态分布表的，但是，在计算机时代，通常都不使用正态分布表了，在R中， pnorm() 这个函数完全可以取代正态分布表了，现在我们使用一个Z-scores的向量来计算一下相应的累积概率，如下所示：
以上就是标准正态分布的 累积分布函数(CDF,Cumulative Distribution Function) 曲线。
简单来说， qnorm 是正态分布 累积分布函数(CDF,Cumulative Distribution Function） 的反函数，也就是说它可以视为 pnorm 的反函数，这里的 q 指的是quantile，即分位数。
使用 qnorm 这个函数可以回答这个问题：正态分布中的第p个分位数的Z-score是多少？
现在我们来计算一下，在正态分布分布中，第50百分位数的Z-score是多少，如下所示：
再来看一个案例：在正态分布中，第96个百分位的Z-score是多少，如下所示：
再来看一个案例：在正态分布中，第99个百分位的Z-score是多少，如下所示：
再来看一下 pnorm() 这个函数，如下所示：
从上面我们可以看到， pnorm 这个函数的功能是，我们知道某个Z-score是多少，它位于哪个分位数上。
接着我们进一步举例来说明一下 qnorm 和 pnorm 的具体功能，如下所示：
现在进行绘图，如下所示：
rnomr() 函数的功能用于生成一组符合正态分布的随机数，在学习各种统计学方法时， rnorm 这个函数应该是最常用的，它的参数有 n , mean ， sd ，其中n表示生成的随机数，mean与sd分别表示正态分布的均值与标准差，现在举个例子，如下所示：
现在我们绘制一下上面的几个向量的直方图，看一下它们的均值是否在70附近，如下所示：
在R语言中，生成不同分布的各种类型的函数都是以d,p,q,r开头的，使用原理跟上面的正态分布都一样。
sample() 函数是一个用于生成随机数的重要的核心函数，如果仅传递一个数值n给它，就会返回一个从1到n的自然数的排列，如果传递是 n:m 就是生成从n到m的随机数，如是是 7,5 ，则会生成5个小于7的随机数，如下所示：
从上面的结果可以看出来，这些数字都是不同的，也就是说，sample函数默认情况下是不重复抽样，每个值只出现一次，如果允许有重复抽样，需要添加参数 replace = TRUE ，如下所示：
sample函数通常会从某些向量中随机挑一些参数，如下所示：
也可以挑日期，如下所示：
上述分布函数前面加上r，p、q、d就可以表示相应的目的：
