Lesson 13: 随机模拟

本课总结

本课学习了一些生成随机数的函数。

  1. sample() 函数,replace 参数默认为 FALSE,意思是不放回取样,如果设为 TRUE,意思则为放回取样。
  2. LETTERS 是 R 中一个预定义的变量,包含所有 26 个英文字母的向量。
  3. rbinom() 是二项式随机行数,用法例 rbinom(1, size = 100, prob = 0.7)
  4. rnorm(10, mean = 100, sd = 25) 正态分布,默认平均值为 0,标准差为 1。
  5. rpois(5, lambda = 10) 泊松分布,生成 5 数字,符合均值为 10 的泊松分布。
  6. my_pois <- replicate(100, rpois(5, 10)),上述操作执行 100 次,生成一个矩阵。
  7. cm <- colMeans(my_pois) ,计算上述矩阵的列平均值。
  8. hist(cm) 查看其分布。

本课详细内容

使用像 R 这样的统计编程语言的一大优点是它有大量的工具来模拟随机数。

本课程假设您熟悉一些常见的概率分布,但这些课题只会在随机数生成方面进行讨论。即使你之前没有这些概念的经验,你也应该能够完成课程并理解主要思想。

我们用来生成随机数的第一个函数是 sample()。使用 ?sample 来调出文档。

让我们来模拟滚动四个六边形的骰子: sample(1:6, 4, replace = TRUE)

> sample(1:6, 4, replace = TRUE)
[1] 1 2 3 4

现在重复该命令,看看结果有何不同。(结果完全相同的概率是 (1/6)的 4 次方,这非常小!)

> sample(1:6, 4, replace = TRUE)
[1] 6 4 4 4

sample(1:6, 4, replace = TRUE) 指示 R 随机选择 1 到 6 之间的四个数字,并可以重复(有放回抽样)。replace 参数的意思是每个数字都是可以重复出现的,因此相同的数字可以出现不止一次,这就是我们想要的,因为您在一个骰子上滚动的结果不会影响您在其他骰子上滚动的结果。

现在在 1 到 20 之间采样 10 个数字,数字不可重复出现。要进行不重复出现的抽样,只需去掉 replace 参数。

> sample(1:20, 10)
 [1] 10 18  1 16  3  5  4 17  7  8

由于最后一个命令在没有 replace 参数的情况下采样,所以输出结果中不会出现重复的数字。

LETTERS 是 R 中一个预定义的变量,包含所有 26 个英文字母的向量。现在就来看看吧。

> LETTERS
 [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L" "M" "N" "O" "P" "Q" "R" "S" "T" "U" "V" "W" "X" "Y" "Z"

还可以使用 sample() 函数对向量的元素进行置换或重新排列。例如,尝试 sample(LETTERS) 来排列英语字母表中的所有的 26 个字母。

> sample(LETTERS)
 [1] "Z" "Y" "X" "D" "O" "S" "V" "E" "I" "W" "U" "G" "B" "C" "L" "K" "R" "F" "T" "M" "Q" "J" "H" "P" "N" "A"

这和从 LETTERS 中抽取 26 个样本是一样的,不可重复出现。当没有指定 sample()size 参数时,R 取一个与您采样的向量大小相同的样本你。

现在,假设我们想要模拟抛 100 次不均匀的两面正面硬币。这个特定的硬币有 0.3 的概率落反面,和 0.7 的概率落正面。

设 0 代表反面,1 代表正面。使用 sample() 从向量 c(0, 1) 中绘制大小为 100 的样本,并可重复出现。因为这枚硬币是不公平的,我们需要第四个参数 prob = c(0.3, 0.7) 来为 0 和 1 赋予特定的概率。将结果分配给一个名为 flips 的新变量,并查看其内容。

> flips <- sample(c(0, 1), 100, replace = TRUE, prob = c(0.3, 0.7))

> flips
  [1] 1 1 1 1 1 1 1 1 1 0 1 1 0 0 0 1 1 1 0 0 1 0 0 1 1 1 1 1 0 0 1 1 1 1 0 0 1 1 1 1 0 0 1 1 0 1 0 0 1 1 1 1 0 1 1 1 0 0 1 1 1
 [62] 0 1 1 1 0 1 1 0 1 0 1 1 1 1 1 1 0 1 0 1 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 0 0 0 1

由于我们将每次抛硬币出现正面的概率设为 0.7,所以我们预计大约有 70 次抛硬币的结果都是 1。使用 sum() 函数计算抛掷中包含的 1 的实际数目。

> sum(flips)
[1] 69

抛硬币是一个二进制结果(0 或 1),我们执行 100 次独立实验(抛硬币),所以我们可以使用 rbinom() 来模拟一个二项随机变量。使用 ?rbinom 调出文档。

R 中的每个概率分布都有一个对应的 r*** 函数(意思是 随机 random),一个 d*** 函数(密度 density),一个 p*** 函数(概率 probability)以及一个 q*** 函数(分位数 quantile)。我们最感兴趣的是本课中 r*** 函数,但我鼓励您自己取探索其他函数。

二项随机变量表示在给定数目的独立试验(抛硬币)中成功(正面)的数目。因此,我们可以使用 rbinom(1, size = 100, prob = 0.7) 生成一个随机变量,它表示非均匀硬币 1001 次抛投中正面的数量。注意,您只需指定成功(正面)的概率,而不需指定失败(反面)的概率。现在试一试。

> rbinom(1, size = 100, pro = 0.7)
[1] 78

同样地,如果我们像看到所以的 0 和 1,我们可以请求 100 个观测值,每个大小为 1,成功的概率为 0.7。尝试一下,将结果分配给一个名为 flips2 的新变量。

> flips2 <- rbinom(size = 1, n = 100, prob = 0.7)

> flips2
  [1] 1 0 0 1 1 0 0 1 1 1 1 0 0 1 1 1 1 0 0 0 1 1 0 1 1 0 0 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 0
 [62] 1 1 1 1 0 0 1 1 1 1 1 0 1 0 0 1 1 0 0 1 1 1 1 1 0 0 1 1 1 0 1 1 0 0 1 0 1 1 1

现在使用 sum() 函数来计算数字 1 出现的次数。它应该很接近 70!

> sum(flips2)
[1] 71

rbinom() 类似,我们可以使用 R 来模拟许多其他概率分布中的随机数。现在打开 rnorm() 的文档。

这是一个均值为 0,标准差为 1 的正态分布。正如您在文档的 Usage 一节中看到的,rnorm()meansd 参数默认值分别是 0 和 1。因此,rnorm(10) 将从一个标准正态分布中产生 10 个随机数。试试吧。

> rnorm(10)
 [1] -1.128544239 -0.582242837  0.006230103 -0.621767254  0.819396336 -2.043350180 -2.085213991  1.865582828 -0.817865961
[10] -1.482506499

现在再做一次同样的事情,只不过平均值是 100,标准差是 25。

> rnorm(10, mean = 100, sd = 25)
 [1]  80.26535 112.38797  97.91502  78.87868  87.42145 122.99880  86.33341 178.35885 134.10795  93.46554

最后,如果我们想要模拟 100 组的随机数,每组包含 5 个由均值为 10 的泊松分布产生的值,我们从一组 5 个数字开始,然后我将向您展示如何以一种方便而紧凑的方式重复 100 次操作。

从均值为 10 的泊松分布中生成 5 个随机值。如果需要帮助,请参阅 rpois() 的文档。

> rpois(5, lambda = 10)
[1]  5 11  7  9 10

现在使用 replicate(100, rpois(5, 10)) 来执行这个操作 100 次。将结果存储再一个名为 my_pois 的新变量里,并查看其内容。

> my_pois <- replicate(100, rpois(5, 10))
> my_pois
     [, 1] [, 2] [, 3] [, 4] [, 5] [, 6] [, 7] [, 8] [, 9] [, 10] [, 11] [, 12] [, 13] [, 14] [, 15] [, 16] [, 17] [, 18] [, 19] [, 20] [, 21]
[1,]   15   12    9    7    8    8    7    7    6    17    16     8    12     8     9    11    11     6     7    12     6
[2,]    8    7    6    9   10    4   13   10    5     7    10     6    14     9    15     7     6    14    12     6     6
[3,]   10    6    6   10   11   10    7    8    9    15     8     7    11    13     3     8     5    12    10     9     3
[4,]   13   11   10    7    8   10   10    7    6     9    13    17    11     4     5    10     9    12     8     9     8
[5,]    9   10   10   11   15    7    2   10   13     8     6     6     7    13    13     4     9    13     9    12    12
     [, 22] [, 23] [, 24] [, 25] [, 26] [, 27] [, 28] [, 29] [, 30] [, 31] [, 32] [, 33] [, 34] [, 35] [, 36] [, 37] [, 38] [, 39] [, 40] [, 41]
[1,]    11     4     6    10     9     5     9     9     7    14    16     6    12     6     9     8    10    10    11     7
[2,]     9    10    13     8     8    19     9     4    17     9    12     8     7     8     9     7     6    11     7    10
[3,]    12    12    13    12     7    11    12     9    10    12    16     7     6    11    10    10     8    11    10    14
[4,]    16    11     9    12    18     8    10     8     9     6     8     9     7    11    13    11     6     8    12    16
[5,]     9    11    12    14    10    12    15    11     9    10    11     5    13     9    13     9     9     9     7     9
     [, 42] [, 43] [, 44] [, 45] [, 46] [, 47] [, 48] [, 49] [, 50] [, 51] [, 52] [, 53] [, 54] [, 55] [, 56] [, 57] [, 58] [, 59] [, 60] [, 61]
[1,]    10     9     9     8     8    13     7     5    11     8    12    10     9    10    18    10     8     7    12    11
[2,]     9    16    13     8     9    16     3    10    11    12    12     9     9    11     8    10    15    12     9    13
[3,]    10     3    15    11    15    11    15    19     8    10    11    13     8     8    10    11     8     6    11    10
[4,]    10     9     9     9     9    17    10     7     9     9    14    14     6     9     7    10    12    12    10     5
[5,]    13    10     8     8    10     8    12    12    11     5     9     9     8    12     9     8     7    12    10    10
     [, 62] [, 63] [, 64] [, 65] [, 66] [, 67] [, 68] [, 69] [, 70] [, 71] [, 72] [, 73] [, 74] [, 75] [, 76] [, 77] [, 78] [, 79] [, 80] [, 81]
[1,]     8    15    10    10    11     6    11     7     7     8    11    13    10     6     4    13     6     7    12     7
[2,]    15     5    16     5     6    11     8     6    15    11    11    11     4     4    11     9     6     5    13    12
[3,]     7     6     6     8     7    11    11     5     4    12     7    10     5     8    19    12    16     4    11     7
[4,]     7    12     7    10     3     9     8    11     6     5     9    13    12     9    11    11     8     7     7    10
[5,]    10     9    24    11     4     6     8    11     8    15     9    14     8    12     8    11     9    14     9     9
     [, 82] [, 83] [, 84] [, 85] [, 86] [, 87] [, 88] [, 89] [, 90] [, 91] [, 92] [, 93] [, 94] [, 95] [, 96] [, 97] [, 98] [, 99] [, 100]
[1,]     7    10    13     5    16    12    13    11    14    12     7    15     8    10     3     5     9    12     12
[2,]    11     7    15    11    13    12    14    11     4    16     7     8    14     8     4     9     9    11     14
[3,]    13     8    12     9    10    12     7    13     8    10     5     7     8    10    10    15    13     9      8
[4,]     6    11    13    11     5     8     6    13     9     9     8    13    19    14    10    11    11     7      7
[5,]     8     9     6     8     9    13     9     8     6    11    10     8    11     6     7    10    11     6     11

replicate() 函数创建了一个矩阵,该矩阵的每一列均包含 5 个随机值,这些随机值是根据 Poisson 分布生成的,且均值为 10。现在,我们可以使用 colMeans() 函数再 my_pois 中找到每一列的平均值。将结果存储再一个名为 cm 的变量中。

> cm <- colMeans(my_pois)

> cm
  [1] 11.0  9.2  8.2  8.8 10.4  7.8  7.8  8.4  7.8 11.2 10.6  8.8 11.0  9.4  9.0  8.0  8.0 11.4  9.2  9.6  7.0 11.4  9.6 10.6
 [25] 11.2 10.4 11.0 11.0  8.2 10.4 10.2 12.6  7.0  9.0  9.0 10.8  9.0  7.8  9.8  9.4 11.2 10.4  9.4 10.8  8.8 10.2 13.0  9.4
 [49] 10.6 10.0  8.8 11.6 11.0  8.0 10.0 10.4  9.8 10.0  9.8 10.4  9.8  9.4  9.4 12.6  8.8  6.2  8.6  9.2  8.0  8.0 10.2  9.4
 [73] 12.2  7.8  7.8 10.6 11.2  9.0  7.4 10.4  9.0  9.0  9.0 11.8  8.8 10.6 11.4  9.8 11.2  8.2 11.6  7.4 10.2 12.0  9.6  6.8
 [97] 10.0 10.6  9.0 10.4

让我们通过使用 hist(cm) 绘制直方图来看一下列的分布。

看起来我们的列均值几乎是正态分布的,对吧这就是中心极限定理,不过这是以后的内容了。

所有的标准概率分布都内置于 R 中,包括指数分布(rexp()),卡方分布(rchisq()), 伽玛分布(rgamma()),…你看到规律了吧。

随机数实际上是一个独立的领域,我们只是触及了它的表面。我鼓励您自己进一步探索这些功能和其他功能。

Copyright © xresearcher.com 2020 all right reserved,powered by Gitbook该文件修订时间: 2020-03-31 15:56:20

results matching ""

    No results matching ""