R官方入门教程（8）概率分布_无法精確計算带连结的p值-CSDN博客

本文链接：https://blog.csdn.net/2301_76324977/article/details/130416649

R是一门著名的可用于数据和统计分析的程序语言，本文翻译自R软件官方文档教程《An Introduction to R》，仅供学习和参考。

8 概率分布

8.1 R作为概率分布表

R 的一个方便用途就是能够提供一套全面的概率分布表。R提供了很多函数来计算概率分布函数 $F (x) = P (X \leq x)$ 、概率密度函数 $f (x)$ 和概率分位数函数（给定 $q$ ，最小的 $x$ 使得 $P (X \leq x) > q$ ）的值，并能够根据分布进行数值模拟。

分布名	R名	可选参数
$B e t a$ 分布	`beta`	`shape1`,`shape2`,`ncp`
二项分布	`binom`	`size`,`prob`
柯西分布	`cauchy`	`location`,`scale`
$\chi^2$ 分布	`chisq`	`df`,`ncp`
指数分布	`exp`	`rate`
$F$ 分布	`f`	`df1`,`df2`,`ncp`
$G amma$ 分布	`gamma`	`shape`,`scale`
几何分布	`geom`	`prob`
超几何分布	`hyper`	`m`,`n`,`k`
对数正态分布	`lnorm`	`meanlog`,`sdlog`
$L o g i s t i c$ 分布	`logis`	`location`,`scale`
负二项分布	`nbinom`	`size`,`prob`
正态分布	`norm`	`mean`,`std`
泊松分布	`pois`	`lambda`
有符号等级分布	`signrank`	`n`
$t$ 分布	`t`	`df`,`ncp`
均匀分布	`unif`	`min`,`max`
威布尔分布	`weibull`	`shape`,`scale`
威尔科克森分布	`wilcox`	`m`,`n`

在R名前加上“d”表示概率密度函数，“p”表示概率分布函数，“q”表示分位数函数，“r”表示随机模拟数据（带有随机偏差）。dxxx的第一个参数是x。pxxx的第一个参数是q。qxxx的第一个参数是p。rxxx的第一个参数是n(除了rhyper,rsignrank,rwilcox,它们是nn)。通常，非中心参数ncp都可用。

pxxx() 和 qxxx() 函数都有对数相关参数 lower.tail 和 log.p ，而 dxxx() 函数有 log。这允许通过-pxxx(t, ..., lower.tail = FALSE, log.p = TRUE) 直接获取累积风险函数值 $H (t)$ ,或通过 dxxx(..., log = TRUE)获取更准确的对数似然函数。
$H (t) = - l o g (1 - F (t))$
此外，函数ptukey()和qtukey()用于正态分布样本的学生化范围的分布，函数dmultinom()和rmultinom()用于多项式分布。更多的分布函数可在各种包中获得，尤其是SuppDists(https://CRAN.R-project.org/package=SuppDists).

下面是一些小例子：

> 2*pt(-2.43, df = 13) # 自由度为13的t分布的双侧p值
[1] 0.0303309
> qf(0.01, 2, 7, lower.tail = FALSE)	# F(2,7)的上1%分位数
[1] 9.546578

8.2 检查一组数据的分布

给定一组（单变量）数据，我们可以通过多种方式检查其分布。最简单的是检查其数字特征。 summary 和 fivenum 给出了两个略有不同的数字特征总结，并用stem()函数显示其数字分布特征（茎叶图）。

> attach(faithful)	# 将R中内置数据集"faithful"加载到工作区
> summary(eruptions)	# 对eruptions进行统计摘要
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.600   2.163   4.000   3.488   4.454   5.100 
> fivenum(eruptions)	# 最小值、第一四分位数、中位数、均值、第三四分位数和最大值
[1] 1.6000 2.1585 4.0000 4.4585 5.1000
> stem(eruptions) # 茎叶图分析

  The decimal point is 1 digit(s) to the left of the |

  16 | 070355555588
  18 | 000022233333335577777777888822335777888
  20 | 00002223378800035778
  22 | 0002335578023578
  24 | 00228
  26 | 23
  28 | 080
  30 | 7
  32 | 2337
  34 | 250077
  36 | 0000823577
  38 | 2333335582225577
  40 | 0000003357788888002233555577778
  42 | 03335555778800233333555577778
  44 | 02222335557780000000023333357778888
  46 | 0000233357700000023578
  48 | 00000022335800333
  50 | 0370

茎叶图很像直方图。R 有一个函数 hist()可以用来来绘制直方图。

> hist(eruptions)

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-xJtMAKqs-1682606056324)(D:\documents\R\notes\assets\image-20230427214303336.png)]$

这样绘制的图横轴划分太粗而且不够具体。我们可以添加一些细节:

> hist(eruptions, seq(1.6, 5.2, 0.2), prob=TRUE)	# 左边界1.6,右边界5.2,间隔0.2
> lines(density(eruptions, bw=0.1))	# 概率密度估计曲线
> rug(eruptions) # 添加边际毛毯

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-s4Kf7BuR-1682606056326)(D:\documents\R\notes\assets\image-20230427214815672.png)]$

我们可以使用函数 ecdf 绘制经验累积分布函数。

plot(ecdf(eruptions), do.points=FALSE, verticals=TRUE)

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-N6VEc2aU-1682606056327)(D:\documents\R\notes\assets\image-20230427215005107.png)]$

这种分布显然不是任何一种标准的分布。切掉一部分数据会不会好一些呢，比如说只留下超过 3 分钟的喷发？让我们试着拟合正态分布。

> long <- eruptions[eruptions > 3]
> plot(ecdf(long), do.points=FALSE, verticals=TRUE)
> x <- seq(3, 5.4, 0.01)
> lines(x, pnorm(x, mean=mean(long), sd=sqrt(var(long))), lty=3)

分位数 - 分位数 (Q-Q) 图可以帮助我们更仔细地验证这一点。

> par(pty="s")	# 方形绘图区
> qqnorm(long) # 绘制与正态分布进行比较的Q-Q图
> qqline(long) # 绘制Q-Q线

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-j4OV7um4-1682606056328)(D:\documents\R\notes\assets\image-20230427220019520.png)]$

上图显示出数据与正态分布具有一定的拟合度，但右尾比正态分布预期的要短。

最后，我们可能想要对数据是否符合正态性进行更正式的检验。 R 提供 Shapiro-Wilk 检验

> shapiro.test(long)

	Shapiro-Wilk normality test

data:  long
W = 0.97934, p-value = 0.01052

和 Kolmogorov-Smirnov 检验

> ks.test(long, "pnorm", mean = mean(long), sd = sqrt(var(long)))

	Asymptotic one-sample Kolmogorov-Smirnov test

data:  long
D = 0.066133, p-value = 0.4284
alternative hypothesis: two-sided

8.3 单样本和双样本检验

到目前为止，我们已经将单个样本与正态分布进行了比较。一个更常见的操作是比较两个样本的各个方面。请注意，在 R 中，所有“经典”测试（包括下面使用的测试）都在通常会被加载的包stat中。

考虑以下从同一实验中不同观测方式得到的两组数据:

Method A: 79.98 80.04 80.02 80.04 80.03 80.03 80.04 79.97
80.05 80.03 80.02 80.00 80.02
Method B: 80.02 79.94 79.98 79.97 79.97 80.03 79.95 79.97

我们对其绘制箱型图以进行比较。

> A <- scan()	# 输入数据
1: 79.98 80.04 80.02 80.04 80.03 80.03 80.04 79.97
9: 
Read 8 items
> B <- scan()
1: 80.02 79.94 79.98 79.97 79.97 80.03 79.95 79.97
9: 
Read 8 items
> boxplot(A, B) # 绘制箱型图

在这里插入图片描述

箱型图表明第一组往往比第二组给出更高的结果。

为了测试两个示例的均值是否相等，我们可以使用未配对的 t 检验：

> t.test(A, B)

	Welch Two Sample t-test

data:  A and B
t = 2.6908, df = 13.823, p-value = 0.01773
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.008078416 0.071921584
sample estimates:
mean of x mean of y 
 80.01875  79.97875

假设数据有一定正态性，确实表现出了显著性。

默认情况下，R函数不假设两个样本的方差相等。我们可以使用F检验来检验方差的相等性，前提是这两个样本来自非异常数据。

> var.test(A, B)

	F test to compare two variances

data:  A and B
F = 0.79673, num df = 7, denom df = 7, p-value = 0.772
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.159509 3.979610
sample estimates:
ratio of variances 
         0.7967332

这并没有显示出显著性，因此我们可以假设方差相等，进行经典的t检验。

> t.test(A, B, var.equal=TRUE)

	Two Sample t-test

data:  A and B
t = 2.6908, df = 14, p-value = 0.01757
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 0.008116742 0.071883258
sample estimates:
mean of x mean of y 
 80.01875  79.97875

所有这些检验都假设两个样本是正态的。而双样本 Wilcoxon（或 MannWhitney）检验仅假设数据分布连续。

> wilcox.test(A, B)
Warning in wilcox.test.default(A, B) : 无法精確計算带连结的p值

	Wilcoxon rank sum test with continuity correction

data:  A and B
W = 54.5, p-value = 0.01902
alternative hypothesis: true location shift is not equal to 0