第9章 R 软件简介及例题

第九章 R 软件简介及例题

§9.1 R 软件简介

R 软件是 S 语言的一种扩展的实现，是一个优秀的统计计算的软件。S 语言最初是 Bell 实验室提出的、主要用于统计分析的交互式语言，S-plus 就是 S 语言商业版本的一个实现，而 R 软件是 S 语言自由软件的一个实现。R 在 GNU 协议 General Public Licence 下免费发行。通常也将 R 所实现的 GNU 协议下的 S 语言扩展称为 R 语言。R 软件最初由 Auckland 大学统计系的 R. Gentleman 和 R. Ihaka 编制，目前这个项目由 R 软件开发核心组维护，其官方网站是 http://www.r-project.org，我们可以通过这个网站了解 R 软件开发的最新进展和 R 软件的最新版本。在 R 软件的官方网站上提供了软件 R 针对 Unix/Linux 和 Windows 等不同操作系统的版本以及大量丰富的扩展包供使用者下载。

R 软件是一个集数据操作、数据分析和图形呈现为一体的免费软件系统，在科研和教育系统得到了广泛认同。一般的软件往往会直接输出分析的结果，而 R 软件则将这些结果都存放在一个个“对象”里面，我们可以有选择的从对象中只抽出我们感兴趣的部分。

读者可以从 R 软件官方网站 http://www.r-project.org 或 CRAN 镜像站点按照所用计算机的操作系统去选择下载相应的 R 软件。目前 R 软件命名形式已统一为 R-x.x.x-win.exe，下载后直接运行，并选择安装目录及选装所需内容即可。安装完毕后，运行 R 软件，你会进入到 R 软件的控制窗口，看到 R 软件就绪等待输入命令的提示符“>”。

R 软件是一种解释型语言，输入的指令可以直接被执行。在提示符后你可以以交互式的行命令方式一个个的输入你的指令，也可以创建一个脚本文件并以批方式运行所写的脚本文件。

许多扩展的、新的开发包可以在 R 软件的菜单界面中 “程序包” 选单中选择安装加载. 例如, 当需要使用机器学习中的支撑向量机方法时, 可以选择安装和加载 e1071 包.

9.1.1 基本操作

在使用 R 软件时，为了便于管理，通常会先创建一些自己的工作目录，用于

在使用 R 软件时保存数据文件。可以简单地利用所用操作系统完成工作目录创建。启动 R 软件后，在还不熟悉 R 软件行命令方式时，有些操作可以通过 R 软件的选单完成。例如在“文件”选单的“改变工作目录”子选单去指定先前所建的工作目录，就可以将原光工作目录转到指定目录。

在 R 软件的提示符 “>” 后，可以输入你想要执行的 R 软件指令。例如，我们希望列出当前目录，>getwd() 每个命令输完后回车，就会看到 R 软件执行的结果。当需要退出 R 软件系统时，可以从 “文件” 选单的子选单 “退出” 选择退出。当然也可以使用行命令方式

q()

此时系统会询问是否保存你目前 R 软件工作空间映像. 如选择保存, 所存储的各类数据在以后的 R 软件任务中可以调出继续使用.

初学者注意 R 软件使用中的一些常识.

(1) 帮助. R 软件有相当丰富的联机帮助文档. 随 R 软件安装附带在帮助选单里的已包括了最常用的几个手册. 官方网站 CRAN 上有大量免费的在线文档可供参考.

(2) R 软件中的函数和关键字. 如果我们知道一个函数和关键字的完整拼写, 例如对于函数 getwd(), 可以使用如下一些方式得到 getwd 的详细说明.

或

如果只记得函数或关键字中的一部分字符，可以使用 apropos 查完整拼写。例如：

只记得 getwd 的 "getw",

apropos("getw")

[1] "getwd" "getWindowsHandle" "getWindowsHandles"

[4] "getWindowTitle" "getWinProgressBar"

此时，R 软件列出了包含 getw 的各条命令。其他常用的还有 find，help.search，args 等。

(3) 指令补全. R 软件提示符“>”后输入指令, 可以在顺序输入部分字符后按“Tab”键补全命令, 如果输入部分字符唯一确定了指令则直接补齐, 若不唯一可再按“Tab”, 将列出含有这些字符的各条指令.

(4) R 软件提供对已输入命令记忆机制，通过使用上、下箭头按键可以在命令记录中向前或向后调出输入过的命令.

9.1.2 数据类型和对象

R 软件所创建、操作的实体是对象. 对象可以是变量、数组、字符串、函数以及由这些元素组成的复杂的结构.

R 软件基本的命令由表达式或者赋值语句组成. 当一个表达式作为一条命令给出时, R 软件会对其求值、打印, 但表达式的值并不保存. 一个赋值语句也会对表达式求值, 并且把表达式的值传给一个变量, 但并不会自动打印出来. 命令间由 “;” 来分隔, 或者另起一行. 注释以井号 (#) 开始, 到行末结束, 不会被执行.

当 R 软件运行时，所有变量、数据、函数及结果都以对象 (objects) 的形式存在计算机的活动内存中。在 R 软件中进行的所有操作都是针对存储在活动内存中的对象。一个对象可以通过赋值操作来产生。在 R 软件中给对象赋值有多种形式，可以直接赋一个数值，也可以是一个算式或一个函数的结果。

可以通过用一些运算符 (如算术、逻辑、比较等) 和一些函数 (其本身也是对象) 来对这些对象进行操作。对象的名字必须是以字母或点 “.” 开头，并且如果对象名开头字母为点时接续的下一个字符不能是数字。对象名中的字符可以是字母 (A ~ Z, a ~ z)、数字 (0 ~ 9)、点 (.) 及下划线 (-). R 软件沿用了 UNIX/Linux 和 Mac OS 的习惯，对象名区分大小写，这样 x 和 X 表示不同的对象。R 软件在涉及路径是也沿用符号 “/”，在 Windows 系统使用 R 软件如涉及路径中 “\” 时可以用 “\\” 代替。R 软件函数总是带有圆括号的形式，即使括号内没有内容 (如, getwd()). 如果直接输入函数名而不输入圆括号，R 软件则会自动显示这个函数的一些具体内容。为了在终端输出对象只需输入对象名即可.

R 软件中所有的对象都有两个内在属性 (intrinsic attributes): 类型和长度。类型是对象元素的基本种类，R 软件中用来表示数据的有数值型、字符型、复数型和逻辑型 (FALSE 或 TRUE, 注意大写) 等。R 软件还可以表示其他一些类型，例如函数。对象的类型和长度可以通过函数 mode() 和 length() 得到。R 软件中不可用或缺失数据 (missing value) 用 NA (not available 不可用) 来表示。R 软件用 Inf 和 -Inf 表示 $ \infty $ 和 $ -\infty $，用 NaN 表示不是数字的值。

当对一个对象不熟悉时，可以用函数 str(), attributes() 等看看该对象里面有什么。函数 str(object) 给出指定对象 object 的结构。函数 attributes(object) 给出指定对象 object 当前定义的非内在属性 (non-intrinsic attributes) 的列表。函数 mode(), length(), attributes(), class() 分别用来给出对象的模式、长度、属性、类别。

赋值：R 软件的赋值操作使用的是“< -”、“- >”和“=”.

运算：算术运算 +,-,*,/,\wedge,\%\,\%/\,,\cdots

运算函数 $ \log(), \exp(), \sin(), \cos(), \tan(), \sqrt{\sin(x)} \cos(x) $,

比较运算 $\gt ,\gt ,=,\lt ,\lt =,-=,-=,!=-,\cdots$.

逻辑运算 !, &, &&, |, ||, xor(), $ \cdots $.

例如：

$ \gt (2+3\times4)\sim2/3\times\sin(pi/2) $

#计算 $ (2+3\times4)^2\div3\times\sin\frac{\pi}{2} $. R 软件中常量 $ \pi $ 用 pi 表示

[1] 65.33333

x<-(2+3*4)^2/3 #将 (2+3×4)^2÷3 的值赋给变量 x，结果并不在终端输出

x<-(2+3*4)^2/3;x #先将 (2+3×4)^2÷3 的值赋给变量 x，之后在终端输出 x

[1] 65.33333

#先将 $ (2+3\times4)^{2}\div3 $ 的值赋给变量 x，输出 x. 同上

[1] 65.33333

#先将 $ (2-3\times4)^{2}\div3 $值赋给变量X，之后输出X，输出x.X≠x

[1] 33.33333

[1] 65.33333

X-x>0 #判断 X-x>0 是否为真

[1] FALSE

1+2+3+ #未输完整回车，注意 R 软件的续行提示符号是 ‘+’

+4+5+6+7+8+9+10

[1] 55

注意, 如果一个对象已经存在, 再对该对象赋值则它以前的值将被新值代替. 表 9.1 给出了 R 软件表示数据的对象的类别概览.

我们主要介绍向量、矩阵、数据框、列表

1. 向量及运算

(1) 向量的生成

建立向量常用函数：

(a) seq() (以及“:”) : 用于建立简单规则的数字向量. rep(): 用于建立规则的向量.

(b) c(), scan(): 多用于建立不规则的一般向量.

例如

x<-c(1,2.2,-1,0); x # 用 1,2.2,-1,0 生成 4 维向量 x，输出

[1] 1.0 2.2 -1.0 0.0

x1<-c("English", "mathematics", "Statistics"); x1 #生成字符型向量 x1，输出

[1] "English" "mathematics" "Statistics"

xs<-seq(10,-5,-pi); xs #生成初始为 10，公差为 -π，终值不小于 -5 的向量

xs，输出

[1] 10.000000 6.858407 3.716815 0.575222 -2.566371

xs1<-pi:3;xs1 #生成初始为 -π，公差为 1，终值不大于 3 的向量

xs1，输出

[1] -3.1415927 -2.1415927 -1.1415927 -0.1415927 0.8584073 1.8584073

2.8584073

Y<-c(x,xs);Y #连接 x,y 构建新向量 Y，输出

[1] 1.000000 2.200000 -1.000000 0.000000 10.000000 6.858407 3.716815

0.575222 -2.566371

Y>0 #给出分量 >0 的真假，逻辑向量

[1] TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE FALSE

rep(1:3,c(3,2,1)) #将（1,2,3）的对应分量分别重复 3，2，1 次. 表达式

[1] 1 1 1 2 2 3

rep(c(1,"c"),3) #将‘'1’，‘'c'’重复三次生成字符向量. 表达式

[1] "1" "c" "1" "c" "1" "c"

rep(c(1,"c"),3)->x3 #将‘'1'，‘'c'重复三次生成字符向量，并赋给 x3

mode(x3) #x3 的类型

[1] "character"

length(x3) #x3 的长度

[1] 6

(2) 向量的运算

R 软件中向量间运算由对应分量运算而得.

c(1, 2, 3, 4)+c(10,20,30,40)

[1] 11 22 33 44

c(1, 2, 3, 4)-c(10,20,30,40)

[1] -9 -18 -27 -36

c(1, 2, 3, 4)*c(10,20,30,40)

[1] 10 40 90 160

c(1,2,3,4)/c(10,20,30,40)

[1] 0.1 0.1 0.1 0.1

sin(c(10,20,30,40))

[1] -0.5440211 0.9129453 -0.9880316 0.7451132

当两个向量维数不等、但较高维是较低维倍数时，将重复较低维向量并运算

1/c(1,2,3,4)

[1] 1.0000000 0.5000000 0.3333333 0.2500000

c(1,2,3,4)+c(100,20)

c(1,2,3,4)+c(100,20,1)

[1] 101 22 4 104

警告信息： $ \mathrm{Inc}(1,2,3,4)+\mathrm{c}(100,20,1) $：长的对象长度不是短的对象长度的整倍数.

(x<--5:10)

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

length(x)

[1] 16

x[3] #x 的第 3 个分量

[1] -3

x[-3] # 向量 x 除去其第三分量后生成的向量

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

x[c(3,4,8)] #向量 x 第 3,4,8 分量生成的向量

[1] -3 -2 2

x[-c(3,4,8)] #向量 x 除去其第 3,4,8 分量后生成的向量

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

x[x<0] #向量 x 小于的分量 0 生成的向量

[1] -5 -4 -3 -2 -1

y<-c(x,"a",1);y #注意 R 软件中数值型向量与字符型向量连接形成的是字符型向量

[1] "-5" "-4" "-3" "-2" "-1" "0" "1" "2" "3" "4" "5" "6" "7"

"8" "9" "10" "a" "1"

2. 矩阵的生成及运算

注意，R 软件中对矩阵还有一个直接扩展：数组 (array, 阵列). 由于比较容

易与矩阵类比，我们不专门介绍.

(1) 矩阵的生成

建立矩阵常用函数：array(), matrix(), cbind(), rbind(), diag(), $ \cdots $.

例如：

x<-1:6; x

[1] 1 2 3 4 5 6

X<-array(x,dim=c(2,3)); X #由 x 生成一个 2 行 3 列的矩阵，注意分量形成的顺序

[,1] [,2] [,3]

[1,] 1 3 5

[2,] 2 4 6

Y<-matrix(x,nrow=3,ncol=2); Y #由 x 生成一个 3 行 2 列的矩阵，注意分量形成的顺序 [,1] [,2]

[1, ] 1 4

[2, ] 2 5

[3, ] 3 6

Xr<-matrix(x,nrow=2,ncol=3,byrow=T); Xr #生成2行3列矩阵，但分量形成的按行顺序[,1] [,2] [,3]

[1, ] 1 2 3

[2, ] 4 5 6

Z<-rbind(X,Xr) #Z 是由 X，Xr 按行合并生成的矩阵

Z

[，1] [，2] [，3]

[1，] 1 3 5

[2，] 2 4 6

[3，] 1 2 3

[4，] 4 5 6

(Z1<-cbind(X,Xr)) #Z 是由 X，Xr 按列合并生成的矩阵

[,1] [,2] [,3] [,4] [,5] [,6]

[1, ] 1 3 5 1 2 3

[2, ] 2 4 6 4 5 6

diag(1:2) #生成对角元为1，2的对角矩阵

[,1] [,2]

[1, ] 1 0

[2, ] 0 2

(2) 矩阵的运算

x<-1:4

[1] 1 2 3 4

(Y<-matrix(x,nrow=2,ncol=2))

[,1] [,2]

[1,1] 1 3

[2,1] 2 4

(Yr<-matrix(x,nrow=2,ncol=2,byrow=T))

[,1] [,2]

[1,1] 1 2

[2,1] 3 4

(Y+Yr) #矩阵加法

[,1] [,2]

[1,1] 2 5

[2,1] 5 8

(Y*Yr) #矩阵逐元相乘

[,1] [,2]

[1,1] 1 6

[2,1] 6 16

(Y%*/Yr) #矩阵乘法

[,1] [,2]

[1,1] 10 14

[2,1] 14 20

(t(Y)) #Y的转置

[,1] [,2]

[1,1] 1 2

[2,1] 3 4

Y[1,2] #Y的第1行、第2列对应元

[1] 3

Y[-1,1] #Y去除其第1行后形成的矩阵

[1] 2 4

W<-matrix(c(4,1,3,1),2,2);W

[,1] [,2]

[1,1] 4 3

[2,1] 1 1

det(W) #矩阵W的行列式

[1] 1

eigen(W) #矩阵W的特征值与特征向量

$values$

[1] 4.7912878 0.2087122

$vectors$

[,1] [,2]

[1, ] 0.9669305 -0.6205203

[2, ] 0.2550401 0.7841904

solve(W) #矩阵W的逆矩阵

[,1] [,2]

[1, ] 1 -3

[2, ] -1 4

qq<-svd(W);qq #W的奇异值分解

$d$

[1] 5.1925824 0.1925824

$u$

[,1] [,2]

[1, ] -0.9628599 -0.2700013

[2, ] -0.2700013 0.9628599

$v$

[,1] [,2]

[1, ] -0.7937170 -0.6082872

[2, ] -0.6082872 0.7937170

class(qq)

[1] "list"

3. 与统计计算有关的常用函数

(1) 与向量和矩阵常联合使用的统计函数

max(), min(), which.max(), which.min(), range(), ...

length(), sum(), prod(), ...

mean(), median(), var(), cov(), cor(), std(),

quantile(), summary(), ...

rev(), sort(), order(), rank(), ...

例如

U<-2:5;U

[1] -2 -1 0 1 2 3 4 5

sum(U)

[1] 12

sum(U)/length(U) #计算 U 的均值

mean(U) #U 的均值

[1] 1.5

U[c(3,4)]<-c(10,7);U #U 的第 3,4 分量值用 10,7 替代

[1] -2 -1 10 7 2 3 4 5

order(U) #U 的元由小到大排列时各元在 U 的原始位置

[1] 1 2 5 6 7 8 4 3

sort(U) #将 U 的元排序

[1] -2 -1 2 3 4 5 7 10

var(U) #U 的方差

[1] 15.71429

summary(U) #U 的统计概要：最小、下四分位数、中位数、均值、上四分位数、最大

Min. 1st Qu. Median Mean 3rd Qu. Max.

-2.00 1.25 3.50 3.50 5.50 10.00

rank(U) #U 中元的秩

[1] 1 2 8 7 3 4 5 6

(2) 常与统计函数联用的两个重要的函数: apply(), sweep() 例如

W<-matrix(c(4,1,3,1),2,2)

(apply(W,MARGIN=2,FUN=mean)) #计算矩阵 W 逐列 (MARGIN=2) 的均值

(FUN=mean)

apply(W, 1, mean) #计算矩阵 W 逐行的均值

[1] 3.5 1.0

row.mean<-apply(W, 2, mean) #计算矩阵 W 逐列均值

sweep(W, MARGIN=2, STATS=row.mean, FUN="--") #将 W 按列减去前步算得的列均值

[,1] [,2]

[1,] 1.5 1

[2,] -1.5 -1

(3) 常用的概率分布有关的函数

注意从指定离散分布总体取样的抽样函数 sample ()，以及组合计数会用到的函数 choose()，factorial()，prod() 等的使用方法.

统计分析软件在统计分析时为处理方便准确，并不像我们手工计算那样提供一个数学用表，通过查表处理，而是提供了更为细致的有关概率函数。这些函数包括四类：d(密度函数或分布列)，p(分布函数)，q(分位数函数)和r(随机数生成函数)。例如：dnorm，pnorm，qnorm 和 rnorm 分别表示正态密度函数、正态分布函数。

布函数、正态分位数函数和正态随机数生成函数. 表9.2中各种分布所含参数的意义从参数称谓上可以了解. 注意，ncp是非中心化参数. 例如

set.seed(321) #指定种子，以便结果重现

u1<-rnorm(1000,mean=3,sd=6) #生成1000个来自均值为3、标准差为6的正态随机数

(mean(u1)) #给出均值估计值

[1] 3.234159

(sd(u1)) #给出标准差估计值

[1] 5.946032

又例如

qnorm(0.975) #给出标准正态分布的0.975分位数

[1] 1.959964

1-pchisq(9.27,3) #对自由度为3, $ \chi^{2} $检验统计量值为9.27时， $ \chi^{2} $检验的p值[1]0.02590836

p 值：统计软件处理假设检验问题时往往直接给出的是 p 值。它所代表的是比由观察值计算出的检验统计值更不利于原假设的可能性的大小。因此，如果一个检验问题显著性为 $ \alpha $，而根据观察计算的检验统计量值算得的 p 值小于 $ \alpha $ 时拒绝原假设，否则不拒绝原假设。

4. 区间估计和假设检验常用函数

(1) t.test(x,...). 单正态总体和两正态总体的均值的 t 检验

t.test(x,y=NULL,

alternative=c("two.sided", "less", "greater"),

mu=0, paired=FALSE, var.equal=FALSE,

conf.level=0.95,...)

主要变量:

x: 非空数值向量.

y: 非空数值向量, 可选 (单总体时省略).

alternative: 指定备择假设类型, 可以是双边检验 ("two.sided", 缺省时), 小于 ("less") 和大于 ("greater") 之一.

mu: 指定期望的真值 (单总体) 或期望值的差 (两总体).

paired: 指定是否是成对样本的 t 检验, 逻辑变量. 缺省时是非成对的检验.

var.equal: 指定两总体的方差是否相等. 缺省时是方差不等的检验.

conf.level: 置信水平. 缺省时是 0.95.

(2) var.test(x,...). 单正态总体和两正态总体的方差的 F 检验

var.test(x,y,ratio=1,

alternative=c("two.sided", "less", "greater"),

conf.level=0.95,...)

主要变量:

x, y: 数值向量, 或拟合线性模型对象.

ratio: x 和 y 的总体方差比的原假设, 默认值是 1.

alternative: 指定备择假设类型, 可以是双边检验 ("two.sided", 缺省时), 小于 ("less") 和大于 ("greater") 之一.

conf.level: 置信水平. 缺省时是 0.95.

(3) chisq.test(). 列联表和拟合优度的 $ \chi^{2} $ 检验

chisq.test(x,y=NULL,correct=TRUE,

p=rep(1/length(x),length(x)),rescale.p=FALSE,

simulate.p.value=FALSE,B=2000)

主要变量:

x: 向量或矩阵.

y: 向量; x 是矩阵时忽略此项.

correct: 逻辑变量, 用于指定是否做连续修正, 默认是 "TRUE".

p: 与 x 有相同维数的概率向量. 默认是均匀分布.

rescale.p: 逻辑变量，指定是否对 p 做归一化，默认是 "FALSE".

simulate.p.value: 逻辑变量, 指定是否用 Monte Carlo 模拟计算 p 值, 默认是 "FALSE".

B: 当 simulate.p.value=TRUE 时, 指定模拟的次数.

5. 列表

R 软件中的列表是由对象的有序集合构成的对象. 列表中包含的对象也称为它的分量 (components). 一个列表的分量可以是不同的模式或类型, 例如一个列表可以同时包括数值向量、逻辑向量、矩阵、复向量、字符数组和函数等等.

例如：

Lst<-list(name="Fred", wife="Mary", no.children=3, child.ages=c(4,7,9))

#创建一个列表

Lst

$name$

[1] "Fred"

$wife$

[1] "Mary"

$no.children$

[1] 3

$child.ages$

[1] 4 7 9

Lst[4] #列表的第4个对象构成的子列表，包括分量对象名

$child.ages$

[1] 4 7 9

Lst[[4]] #列表 Lst 的第4个分量对象

[1] 4 7 9

Lst[[4][2]] #列表 Lst 的第4个分量的第2分量

[1] 7

6. 数据框

数据框可看做是一个由不同模式和属性的列构成的“矩阵”。它常以类似矩阵的形式出现，但可以由不同模式和属性列组成。

数据框的列表对象有下面一些限制条件：

(a) 分量必须是向量 (数值型, 字符型, 逻辑型), 因子, 数值矩阵, 列表或者其他数据框.

(b) 矩阵、列表和数据框成员为新数据框提供了和其列数、成员数、变量数相同个数的变量.

(c) 数值向量、逻辑值、因子保持原有格式，而字符向量会被强制转换成因子并且它的水平就是向量中出现的独立值.

(d) 在数据框中以变量形式出现的向量结构必须长度一致，矩阵结构必须有一样的行数.

例如：

C<-c("1","1","2","3","3","2","3","1")

S<-c(65,57,78,39,87,89,99,92)

G<-c("M","F","F","M","M","F","F","F")

Dsg<-data.frame(class=C,score=S,gender=G); Dsg

Dsg[2,3]

[1] F

Levels: F M

Dsgs<-subset(Dsg,gender="F"); Dsgs

Class score gender

2 1 57 F

3 2 78 F

6 2 89 F

7 3 99 F

8 1 92 F

9.1.3 数据读写

如果我们指定的工作目录是 windows 系统下 d:\R-cx，可以先改变工作目录。例如

getwd() # 查看当前目录

[1]"C:/Documents and Settings/toshiba/My Documents"

setwd("d:\R-ex") #或 setwd("d:/R-ex"), 转到目录 d:\R-ex

getwd()

[1] "d:/R-ex"

假如. 有 Excel 文件在 Excel 中显示如下:

我们可以在 Excel 里的 “文件” 选单 “另存为” 子选单, 例如选择逗号分隔的.csv 格式, 存为 d:\R-ex\ex1.csv. 可以用记事本打开, 看到

为了在 R 软件中读入 d:\R-ex\ex1.csv,

ex1<-read.table("d:/R-ex/ex1.csv",sep="",header=TRUE) #注意 ex1.csv 有标题行

ex1

学号 数学 语文 英语 班级

1 x01 83 89 88 1

2 x02 80 98 90 1

3 x03 95 85 91 2

dim(ex1)

[1] 3 5

attributes(ex1)

$names$

[1] "学号" "数学" "语文" "英语" "班级"

$class$

[1] "data.frame"

$row.names$

[1] 1 2 3

str(ex1)

'data.frame': 3 obs. of 5 variables:

$ 学号: Factor w/ 3 levels "x01","x02","x03": 1 2 3$

$ 数学: int 83 80 95$

$ 语文: int 89 98 85$

$ 英语: int 88 90 91$

$ 班级: int 1 1 2$

ex2<-read.csv("d:/R-ex/ex1.csv",sep="",header=TRUE)

write.csv(ex1,"d:/R-ex/wex1.csv") #将 ex1 以.csv 格式写入 d:/R-ex/wex1.csv

#为保存对象 x 和 y 到一个 RData 文件 xy.RData，可执行

save(x, y, file = "xy.RData") #为加载原先保存的一个 RData 文件 xy.RData，可执行

load('xy.RData')

注意：

(a) 函数 attach() 可以用来连接数据框. R 软件中编辑数据框的个别修改时, 用 fix() 比较方便.

(b) 如果需要保存当前工作空间所有对象，可以在控制台中“文件”选单的“保存工作空间”，指定文件路径和文件名，并“保存”。今后可以从“文件”选单的“加载工作空间”中指定你想加载的文件，这将可以还原指定的原先工作空间。R 软件中也提供了执行这些处理的行命令 save.image()。

9.1.4 简单编程

1. 控制流

R 软件是一个表达式语言，其任何一个语句都可以看成是一个表达式。表达式间以分号分隔或用换行分隔。表达式可以续行，只要前一行不是完整表达式则下一行为上一行的继续。若干个表达式可以放在一起组成一个复合表达式，作为一个表达式使用。组合用花括号“{ }”表示。R 软件也提供了分支、循环等程序控制结构。

(1) 分支语句

分支语句有 if / else 语句、switch 语句.

(a) if / else 语句的格式为

if(cond_1)

statement_1

else if(cond_2)

statement_2

else

statement_3

例如：

x<-c(-2,2,0)

if(any(x<=0)) y<-log(1+x^2) else y<-log(x)

y

[1] 1.609438 1.609438 0.000000

(b) switch 语句

switch 语句是多分支语句，其使用方法是

switch (statement, list)

其中 statement 是表达式，list 是列表，可以用有名定义。如果表达式的返回值在 1 到 length(list)，则返回列表相应位置的值；否则返回“NULL”值。例如：

x<-3

switch(x, 1:3, mean(1:10), sd(2:6))

[1] 1.581139

switch(4, 1:3, mean(1:10), sd(2:6))

NULL

当 list 是有名定义时, statement 等于变量名时, 返回变量名对应的值; 否则返回 “NULL” 值.

(2) 中止语句与空语句

中止语句是 break 语句，break 语句的作用是中止循环，使程序跳到循环以外。空语句是 next 语句，next 语句是继续执行，而不执行某个实质性的内容。

(3) 循环语句

循环语句有 for 循环、while 循环和 repeat 循环语句.

(a) for 循环语句

for 循环的格式为

for (name in expr_1) expr_2

其中 name 是循环变量, expr_1 是一个向量表达式, expr_2 通常是一组表达式.

例如

x<-array(0, dim=c(2,3))

for (i in 1:2){

for(j in 1:3){

x[i,j]<-i+j

}

}

x

[,1] [,2] [,3]

[1,] 2 3 4

[2,] 3 4 5

(b) while 循环语句

while 循环语句的格式为

while (condition) expr

当条件 condition 成立. 则执行表达式 expr.

例如

a<-4

while(a>0){

if(a==2){

cat("before break:", a,"\n")

break

}

a<-a-1

cat(a,"\n")

}

结果显示为

before break: 2

(c) repeat 循环语句

repeat 语句的格式为

repeat expr

repeat 循环依赖 break 语句跳出循环.

例如，计算 1000 以内的 Fibonacci 数.

f<-1; f[2]<=-1; i<-1

repeat{

f[i+2]<=-f[i]+f[i+1]

i<-i+1

if(f[i]+f[i+1]>=1000) break

}

或将条件语句改为 if(f[i]+f[i+1]<1000) next else break, 也有同样的计算结果.

注意，R 软件中循环是相当耗时的，如有可能尽量转化为等价的向量表示.

2. 编写函数

R 软件允许用户自己创建函数. 有许多 R 软件函数存贮为特殊的内部形式，并可以被进一步的调用.

函数定义的格式如下，

name<-function(arg_1, arg_2, ...) expression

expression 是 R 软件中的一组表达式, arg_1, arg_2, … 表示函数的参数. 表达式中, 放在程序最后的信息是函数的返回值, 返回值可以是向量、数组 (矩阵)、列表或数据框.

调用函数的格式为 name(expr_1, expr_2, ...), 并且在任何时调用都是合法的.

在调用自己编写的函数(程序)时，需要将已写好的函数调到内存中. 如果自己编写的程序已存为. R文件，可以使用“运行R脚本文件”选单或执行source()函数运行.

例 9.1.1 写一个利用方差求标准差的函数.

std.dev <- function(x) sqrt(var(x))

计算 t 检验双边 p 值的函数，调用先前所写的求标准差的函数.

t.test.p<-function(x, mu=0)

{

n<-length(x)

t<-sqrt(n)*(mean(x)-mu) / std.dev(x)

2*(1-pt(abs(t), n-1))

}

运行的例子

x<-rnorm(100)

t.test.p(x)

[1] 0.1373749

t.test.p(x,1) #注意，结果说明了什么？

[1] 6.163958e -13

例 9.1.2 计算样本中位数值

s.median<-function(x){

odd.even<-length(x)%2

if(odd.even==0){

(sort(x)[length(x/2)]+sort(x)[1+length(x)/2])/2

}else{

sort(x)[ceiling(length(x)/2)]

}

}

运行的例子

x<-rnorm(100)

s.median(x)

[1] 1.378735

9.1.5 绘图

R 软件中绘图语句分为三个基本类别

(b) 低层 (low-level) 函数：绘图函数在已有图形上添加更多信息，例如额外的点、线和标签.

(a) 高层 (high-level) 函数：绘图函数在图形设备上创建一个新图形，通常包括坐标轴，标签，标题等等.

(c) 交互（Interactive）函数：图形函数允许用户通过鼠标一类的指点设备向已有图形交互的增加信息，或者从中释放信息.

此外，R 软件包含一系列可以用来定制图形的图形参数.

作图时常用函数

windows(): 打开一个绘图窗.

par(mfrow=c(m,n)): 一页绘制 m × n 张图.

1. 高层绘图命令

(1) 高层绘图函数，由函数参数提供数据生成一幅完整的图形。其中适当的坐标轴、标签和标题都可自动生成。

(a) plot()

plot(x,y): 如果 x, y 是向量, plot(x,y) 生成一幅 y 对 x 的散点图. 如果 x 是一个因子对象, y 是一个数值型向量, 则生成 y 对应于 x 各个水平的箱线图.

plot(x): x 是包含两个变量的列表或一个两列的矩阵, 这个命令生成 y 对 x 的散点图; 如果 x 是一个时间序列, 这个命令生成一个时序图; 如果 x 是一个数值型向量, 则生成一个向量值对它下标的散点图; 而如果 x 是一个复向量, 则生成一个向量中元素的虚部对实部的散点图. 如果 x 是一个因子对象, 生成 x 的条形图.

(b) pairs(X), coplot()

若 X 是一个数值矩阵或数据框，则 pairs(X) 生成一个配对的散点图矩阵，矩阵由 X 中的每列的列变量对其他各列列变量的散点图组成，得到的矩阵中每个散点图行、列长度都是固定的。当问题涉及三、四个变量时，使用 coplot() 更好些。如果 a 和 b 是数值向量，c 是数值向量或因子对象（全都是相同长度的），则命令 coplot(a～b|c) 对应 c 的某些给定值生成数个 a 对 b 的散点图。

(c) 其他高层图形函数生成不同类型的图形. 例如:

qqnorm(x): 正态分位数 - 分位数图:

qqline(x): 绘一条过正态分布分位数和数据分位数的第一和第三四分位点的直线.

qqplot(x,y): y 对 x 的分位数 - 分位数图.

hist(x): 直方图.

还有许多, 如 dotchart(x, ...), image(x, y, z, ...), contour(x, y, z, ...), persp(x, y, z, ...), ...

(2) 高层绘图函数的参数

高级图形函数可以使用一系列的参数. 例如:

add=TRUE：强制函数按照低层图形函数的方式操作，将图形置于当前图形上.

axes=FALSE：暂时禁止坐标轴的生成，以便使用 axis() 函数添加你自己定制的坐标轴． 默认情况是 axes=TRUE，即包含坐标轴．

log="x"、log="y"、log="xy"：令 x，y 或者两者全都对数化.

type=参数：type 参数控制所生成图形的类型．主要有

type="p": 绘点（默认值）.

type="l": 绘线.

type="b": 绘点并用线连接.

type="o": 将点绘在线上.

type="h": 绘制从点到横轴的垂线.

type="s": 阶梯式图.

type="n": 不画任何点、线，但绘出坐标轴（默认情况），而且要根据数据绘出坐标系，为给后续的低级图形函数创建图形作基础.

xlab=string，ylab=string：x 轴或 y 轴的标签．使用这些参数来改变默认的标签，通常的默认值是调用高层绘图函数时所使用对象的名称．

main=string: 图表标题，位于图形的顶部，大字体显示.

sub=string: 子标题，位于 x 轴下面，用较小的字体显示.

例如：

par(mfrow=c(1,2))

x<-50:1; y<-x+runif(50, min=10, max=20)

plot(x, y)

plot(x,y,type="1")

运行的结果如图 9.1 所示.

2. 低层绘图命令

有时高层绘图函数并不能很精确生成我们想要的图形. 此时, 我们可以通过低层绘图命令在当前图形上添加信息 (例如, 点、线或文本).

下面列出了一些比较有用的低级绘图函数

points(x, y) 和 lines(x, y) 用于在当前图形上添加点或线. 函数 plot() 的参

数 type 也可以用于这些函数 (默认的是 "p" 代表 points() 和 "l" 代表 lines()). text(x, y, labels, ...) 给定点坐标 x, y, 在该点添加文本. 通常 labels 是一个整数或字符向量, 其中 labels[i] 出现在点 (x[i], y[i]). 默认值是 1 : length(x).

abline(a, b): 在当前图上添加一条斜率为 b，截距为 a 的直线.

abline(h = y): 在图形指定的高度上绘制一条贯穿图形的水平线.

abline(v = x): 在 x 轴的指定位置绘制一条贯穿的垂线.

polygon(x, y, ...): 绘制一个多边形, 其顶点由 (x, y) 指定. 同时还 (可选的) 可以加上阴影线, 如果图形设备允许的话还可以将多边形填充.

legend(x, y, legend, ...): 在当前图形的指定位置添加图例。绘制的字符，线条类型，颜色等等由字符向量 legend 指定。

title(main, sub): 在当前图形的顶部用大字题添加一个标题 main, 在底部用较小的字体添加子标题 sub.

axis(side, ...): 在当前图形的指定边上添加坐标.

低层绘图函数通常都需要一些位置信息（例如，x，y 坐标）来决定在哪里添加新的元素.

数学注释：某些情况下需要在图形中加入数学符号或公式。在 R 软件中可以通过在 text, mtext, axis 或 title 中指定一个表达式来实现。例如，下面的代码绘制了二项概率函数的公式：

例如：

x<-seq(-pi,pi,len=65)

plot(x,sin(x),type="1",ylim=c(-2,3),col=3,lty=2)

points(x,cos(x),pch=3,col=4)

lines(x,tan(x),type="b",lty=1,pch=4,col=6)

legend(-1,3,c("sin", "cos", "tan"), col=c(3,4,6), text.col="green4", lty=c(2,-1,1), pch=c(-1,3,4), bg='gray90')

运行的结果如图 9.2 所示.

§9.2 例题

本节列出的 R 软件程序代码（简称 R 代码）往往只是按容易想到的做法编出的，未经过优化的。读者可以再修改，使其更精炼有效。为节省篇幅，一般不列出运行结果。

例 9.2.1 一副扑克去掉大小鬼，从剩下的这 52 张中任取 5 张，求下列事件的概率：

(1) 同花顺 5 张.

(2) 异花顺 5 张.

(3) 5 张牌中只有两个点数.

(4) 5 张牌中恰有两对.

解 记问题 (1), (2), (3), (4) 对应的事件分别为 $ A_{1}, A_{2}, A_{3}, A_{4} $. 从 52 张牌中任取 5 张, 不计次序有 $ \binom{52}{5} $ 种取法.

(1) 52 张牌中有 4 个花色，每个花色 13 张牌，而每个花色能构成的 5 张同花顺只有 9 种情况。取 5 张构成同花顺有 $ \binom{4}{1} \binom{9}{1} $ 种方法，故

$$ P(A_{1})=\frac{\binom{4}{1}\binom{9}{1}}{\binom{52}{5}}. $$

(2) 先计算 5 张牌构成顺子的数目. 顺子的 5 个点数可由最小点数确定、有 9 种取法, 这 5 个点数每个有 4 个可能的花色, 故可以构成 $ \begin{pmatrix}9\\1\end{pmatrix}4^{5} $ 个顺子. 其中同花顺有 $ \begin{pmatrix}4\\1\end{pmatrix}\begin{pmatrix}9\\1\end{pmatrix} $ 个. 故异花顺有 $ \begin{pmatrix}9\\1\end{pmatrix}\begin{bmatrix}4^{5}-(4)\\1\end{bmatrix} $ 个, 于是

$$ P(A_{2})=\frac{\binom{9}{1}\left[4^{5}-\binom{4}{1}\right]}{\binom{52}{5}} $$

(3)5张牌中只有两个点数,必定是有3张同一个点数、另两张是另外一个点数,考虑到扑克有4个花色,故5张牌中只有两个点数有 $ \left[\binom{13}{1}\binom{4}{3}\right]\left[\binom{12}{1}\binom{4}{2}\right] $种可能.故

$$ P(A_{3})=\frac{\left[\binom{13}{1}\binom{4}{3}\right]\left[\binom{12}{1}\binom{4}{2}\right]}{\binom{52}{5}} $$

(4) 考虑到两个对子有 $ \left[\binom{13}{1}\binom{4}{2}\right]\left[\binom{12}{1}\binom{4}{2}\right]/2 $ 种取法，剩下一个单张有 $ \left[\binom{11}{1}\binom{4}{1}\right] $ 种取法，故

$$ P(A_{4})=\frac{\left[\binom{13}{1}\binom{4}{2}\right]\left[\binom{12}{1}\binom{4}{2}\right]\bigg/2}{\binom{52}{5}}=\frac{2\;A_{13}^{3}\left[\binom{4}{2}\right]^{2}}{\binom{52}{5}}. $$

R 代码：

PA1<-4*9/choose(52,5); PA1 #理论值

set.seed(330311)

x<-rep(1:4,c(13,13,13,13)) #模拟计算 PA1

y<-rep(1:13,4)

n<-300000

S<-cbind(X=x,Y=y)

k1<-numeric(n)

for(i in 1:n)

{T<-sample(1:52,5,replace=FALSE)

u<-S[T,]

k1[i�<-(max(u[,1])-min(u[,1])<1)&(max(u[,2])-min(u[,2])==4))}

list(PA1 模拟 =sum(k1)/n) #模拟

PA2<-9*(4~5-4)/choose(52,5); PA2 #理论值

PA3<-13*12*choose(4,3)*choose(4,2) /choose(52,5); PA3 #理论值

PA4<-prod(13:11)*choose(4,2)^2*2/choose(52,5); PA4 #理论值

例 9.2.2 假设每个人的生日出现在一年 365 天中的任一天是等可能的.

(1) 生日问题: 随机抽取 n 个人组成一群体, 问这个群体至少有两个人生日同一天的概率有多大? 并作出随 n 不同概率取值的曲线.

(2) 给定一年中的特定一天，例如十月一日，为使随机抽取 n 个人组成群体中有某人的生日恰为这一天的概率超过 0.5，n 至少取多大？

（3）随机抽取3人 $ S_{1}, S_{2}, S_{3} $，记事件 $ A=\{S_{1} $与 $ S_{2} $生日是同一天\}, $ B=\{S_{2} $与 $ S_{3} $生日是同一天\}, $ C=\{S_{3} $与 $ S_{1} $生日是同一天\}, $证明A,B,C 两两独立，但不相互独立.$

$$ 解 \quad\Omega=\{(m_{1},m_{2},\cdots,m_{n})|m_{i}\in\{1,2,\cdots,365\},i=1,2,\cdots,n\}. $$

A = {至少有两个人生日相同}.

D = {恰有同学生日是十月一日}.

(1) $P(\overline{A})=\frac{A_{365}^{n}}{365^{n}}=\frac{365\times364\times\cdots\times(365-n+1)}{365^{n}}$, $P(A)=1-\frac{A_{365}^{n}}{365^{n}}$.

(2) $P(D)=1-\left(1-\frac{1}{365}\right)^{n}=1-\left(\frac{364}{365}\right)^{n}\geqslant0.5$, $n\geqslant\frac{\log0.5}{\log\frac{364}{365}}\approx252.652$.

(3) $P(A)=P(B)=P(C)=\frac{\binom{365}{1}}{365^{2}}=\frac{1}{365}$, $P(AB)=P(BC)=P(CA)=\frac{\binom{365}{1}}{365^{2}}=\frac{1}{365^{2}}$, $P(ABC)=P(AB)=\frac{1}{365^{2}}$, $P(A)P(B)P(C)=\frac{1}{365^{3}}\neq P(ABC)$,

可见 A, B, C 两两独立，但不相互独立.

R 代码：

m<-1:100 #m<-c(1,seq(2,100,3))

p.bir<-function(m){c(m,1-prod(365:(365-m+1))/365^m)}

t.bir<-t(supply(m,p.bir))

colname1<-c("m","prob")

dimnames(t.bir)<-list(NULL, colname1)

x<-c(1,seq(5,60,5),seq(65,100,5))

t.bir[x,]

plot(t.bir,

xlab = "群体中人数", ylab = "群体中至少两人同天生日的概率",

#main = "生日问题",

xlim=c(0,100), ylim=c(0,1),

xact="n", yact="n")

axis(1,at= seq(0,100,5), labels=as.character(seq(0,100,5)),las=1)

axis(2,at=seq(0.1,1,0.1), labels=as.character(seq(0.1,1,0.1)),las=1)

lines(t.bir)

abline(h=0.5)

abline(v=23, lty=2) # dashed line

floor(log(0.5)/log(364/365))+1

运行结果中的图如图 9.3 所示.

例 9.2.3 一批产品 10 件中有两件次品、8 件正品. 现每次从中任取一件，且取后不放回，试求下列事件的概率.

(1) 前两次取得正品.

(2) 第二次取得次品.

(3) 若已知第二次取得次品，第一次取得次品.

解 记 $ A_{i} $ 为事件“第 i 次取到次品”，i=1,2.

$$ (1)~P(\overline{A_{1}}~\overline{A_{2}})=P(\overline{A_{1}})P(\overline{A_{2}}|\overline{A_{1}})=\frac{8}{10}\times\frac{7}{9}=\frac{28}{45}. $$

(2) 由全概率公式有

$$ P(A_{2})=P(A_{1})P(A_{2}|A_{1})+P(\overline{A_{1}})P(A_{2}|\overline{A_{1}})=\frac{2}{10}\times\frac{1}{9}+\frac{8}{10}\times\frac{2}{9}=\frac{1}{5}. $$

(3) 由贝叶斯公式有

$$ P(A_{1}|A_{2})=\frac{P(A_{1})P(A_{2}|A_{1})}{P(A_{2})}=\frac{\frac{2}{10}\times\frac{1}{9}}{\frac{1}{5}}=\frac{1}{9}. $$

R 代码：

$$ (P1\lt -8*7/(10*9)) $$

$$ \left(\mathrm{P}2\lt -\left(2*1+8*2\right)/\left(10*9\right)\right) $$

$$ \left(\mathrm{P}3\lt -2*1/\left(10*9\right)/\left(1/5\right)\right) $$

例 9.2.4 甲盒中有 3 个白球、2 个黑球，乙盒中有 3 个白球、1 个黑球，丙盒中有 1 个白球、3 个黑球。现在从甲盒中任取一球放入乙盒，再从乙盒中任取一球放入丙盒，然后从丙盒中任取一球，问最后这个从丙盒中取到的是白球的概率是多少？

解 分别记事件 A, B, C 为从甲、乙、丙取到的是白球. 则

$$ P(A)=\frac{3}{5},\quad P(\overline{A})=\frac{2}{5}, $$

$$ P(B)=P(A)P(B|A)+P(\overline{A})P(B|\overline{A})=\frac{3}{5}\times\frac{3+1}{4+1}+\frac{2}{5}\times\frac{3}{4+1}=\frac{18}{25}, $$

$$ P(C)=P(B)P(C|B)+P(\overline{B})P(C|\overline{B})=\frac{18}{25}\times\frac{2}{5}+\left(1-\frac{18}{25}\right)\times\frac{1}{5}=\frac{43}{125}.\quad\Box $$

R 代码：

$$ \mathrm{Pb}\lt -{\mathrm{Pa}}*4/5+(1-\mathrm{Pa})*3/5 $$

$$ \mathrm{P c}\lt -{\mathrm{P b}}*2/5+({\mathrm{1-P b}})*1/5;\mathrm{P c} $$

例 9.2.5 甲袋中装有 2 个白球、1 个黑球，乙袋中装有 3 个白球。现从甲、乙两袋中任抽一个、交换放入另一袋，这种交换进行了 n 次后，问此时黑球在甲袋、还是在乙袋的可能性大？

解 由于只有一个黑球，若将每次交换球时抽到黑球看作成功，成功率是 p，则 n 次交换可当做是 n 重伯努利实验，事件 A = {n 次交换后黑球在甲袋中}，则 A = {n 次实验中成功了偶数次}。而 “n 重伯努利实验中成功了偶数次” 的概率为

$$ P=\sum_{k=0\bmod2}\binom{n}{k}p^{k}q^{n-k}, $$

其中 q = 1 - p. 但

$$ 1=(p+q)^{n}=\sum_{k=0\bmod2}{\binom{n}{k}}p^{k}q^{n-k}+\sum_{k=1\bmod2}{\binom{n}{k}}p^{k}q^{n-k}, $$

$$ (q-p)^{n}=\sum_{\substack{k=0\mathrm{~m o d~}2}}\binom{n}{k}p^{k}q^{n-k}-\sum_{\substack{k=1\mathrm{~m o d~}2}}\binom{n}{k}p^{k}q^{n-k}, $$

故

$$ P=\frac{1}{2}[1+(p-q)^{n}]=\frac{1}{2}[1+(1-2p)^{n}]. $$

对本题 p = 1/3, 于是 $ P = \frac{1}{2} \left(1 + \frac{1}{3^n}\right) $, 说明黑球在甲袋的可能性大些. □

R 代码：

ex5.p1<-function(p=1/3,n=10){ #理论值

p1<-1/2*(1+(1-2*p)^n);p1

}

ex5.p1(,6)

set.seed(330311)

ex5.p11<-function(p=1/3,n=10,m=1000){ #模拟

x<-numeric(m)

for (j in 1:m)

{x[j]<sum(sample(0:1,n,replace=TRUE))%%2}

p11<-sum(x==0)/m;p11

}

ex5.p11()

例 9.2.6 某厂组装的器件有 70% 可以直接出厂销售，剩下的 30% 需进行调试，这些经调试的器件经再测试其中有 80% 可以出厂销售，另 20% 成为废品。假定该厂生产了 96 个器件，且各器件能否出厂销售相互独立，试求以下事件的概率。

(1) 这批产品都可出厂销售.

(2) 这批产品中至少有两件废品.

解 对该厂的任一器件，记 A：“可直接出厂销售”，B：“可出厂销售”，则

$$ P(B)=P(A)P(B|A)+P(\overline{A})P(B|\overline{A})=0.7\times1+0.3\times0.8=0.94. $$

(1) $ p_1 = \begin{pmatrix} 96 \\ 2 \end{pmatrix} \times 0.94^{94} \times (1 - 0.94)^2 $.

$$ (2)\;p_{2}=1-[0.94^{96}+\binom{96}{1}\times0.94^{95}\times(1-0.94). $$

R 代码：

ex6.p<-function(n,r,s)

{

p<-r*1+(1-r)*s

p1<-choose(n,2)*p^-(n-2)*(1-p)^2

p2<-1-(p^n+n*p^-(n-1)*(1-p))

list(P1=p1,P2=p2)

}

ex6.p(96,0.7,0.8) #n<-96 p<-0.94

例 9.2.7 Monty Hall 问题. 这是一个有奖竞猜问题. 主持人在现场准备一些道具和奖品，有3扇闭合的门，1件大奖汽车和2件小奖山羊. 主持人将这三件奖品随机分派在这三个门后，让参与者猜大奖藏在哪个门后，当参与者选择一个门后，主持人会从参与者选择的门之外的其他两扇门中随意打开一扇没有大奖的门给参与者看，然后允许参与者改换先前的选择或不换，主持人以这时参与者确定的选择为准，将门后的奖品奖给参与者. 问参与者在主持人允许改换时选择“换门”、还是“不换门”得大奖的可能性更大?

解　为简便起见，将参与者首次选择的门记为1号门，其他两扇分别记为2，3号门， $ A_{i}=\{ $大奖在第i号门 $ \} $， $ B_{j}=\{ $参与者首次选择后主持人打开的是第j门 $ \} $，则

$$ P(A_{1})=P(A_{2})=P(A_{3})=\frac{1}{3}, $$

现若主持人打开的是2号门，则

$$ P(B_{2}|A_{1})=\frac{1}{3},\quad P(B_{2}|A_{2})=0,\quad P(B_{2}|A_{3})=1. $$

“不换门”得大奖的概率为

$$ \begin{align*}P(A_{1}|B_{2})&=\frac{P(A_{1})P(B_{2}|A_{1})}{P(A_{1})P(B_{2}|A_{1})+P(A_{2})P(B_{2}|A_{2})+P(A_{3})P(B_{2}|A_{3})}\\&=\frac{\frac{1}{3}\times\frac{1}{2}}{\frac{1}{3}\times\frac{1}{2}+\frac{1}{3}\times0+\frac{1}{3}\times1}=\frac{1}{3},\end{align*} $$

也就是说，此时参与者“换门”选剩下的第3号门中大奖的概率将是 $ 1-1/3=2/3 $。故“换门”中大奖的可能性更大。类似地，参与者首选后主持人打开是3号门时，结论也是选择“换门”中大奖的可能性更大。

R 代码：

set.seed(310311)

example.monty<-function(n=10000){ #模拟 n<-10000

actu<-sample(1:3,n,replace=TRUE)

gues<-sample(1:3,n,replace=TRUE)

equa<-(actu==gues)

pNoSwitch<-sum(equa)/n

notEq<-(actu!=gues)

pSwitch<-sum(notEq)/n

Probs<-c(pNoSwitch,pSwitch)

names(Probs)<-c("P(不换获大奖)","P(换获大奖)")

例 9.2.8 从 1 至 100 这 100 个数中，任取一数。求如下事件的概率。

(1) 取到的数能被4整除.

(2) 取到的数能被6整除.

(3) 取到的数能被4和6整除.

(4) 取到的数能被4或6整除.

(5) 取到的数能被4整除但不能被6整除.

解 记事件 $ A=\{ $取到的数能被4整除 $ \} $， $ B=\{ $取到的数能被6整除 $ \} $.

(1) 由于 $ 100 = 25 \times 4 $，故 1 至 100 中有 25 个数能被 4 整除，于是 $ P(A) = \frac{25}{100} = \frac{1}{4} $.

(2) 由于 $ 100 = 16 \times 6 + 4 $，故 1 至 100 中有 16 个数能被 6 整除，于是 $ P(B) = \frac{16}{100} = \frac{4}{25} $.

(3) $ 4 = 2 \times 2, 6 = 2 \times 3 $，能被4且被6整除，即能被12整除。 $ 100 = 8 \times 12 + 4 $，故 $ P(AB) = 2/25 $.

$$ (4)P(A\cup B)=P(A)+P(B)-P(AB)=1/4+4/25-8/100=33/100. $$

$$ (5)P(A-B)=P(A)-P(AB)=1/4-8/100=17/100. $$

R 代码：

ex8.p1<-function(n=100,a=4,b=6){

p<-numeric(5)

x<-1:n

p[1]<~sum(x%%4==0)/n

p[2]<~sum(x%%6==0)/n

p[3]<~sum((x%%6==0)&(x%%4==0))/n

p[4]<~sum((x%%6==0)|(x%%4==0))/n

p[5]<~sum((x%%4==0)&(x%%6!=0))/n

p

}

ex8.p1()

例 9.2.9 某专业今年计划通过考试招研究生 25 人，确定此次考试总分满分为 500 分，按成绩由高到低安排录取，其中 10 名免费生和 15 名收费生。考试阅卷后该专业公布的考试情况是有 512 人参加了考试，考生平均成绩 280，并划定分数线为 320。现有考生得知自己总分为 333，问该考生有没有可能被录为免

费生？

解 假定考试总分 X 服从正态分布, $ X \sim N(\mu, \sigma^{2}) $.

由设知 $ \mu = 320 $. 近似地有

$$ P(X\geqslant320)\approx\frac{25}{512}=0.04883, $$

由此

$$ 1-\varPhi\left(\frac{320-280}{\sigma}\right)\approx0.04883,\quad\varPhi\left(\frac{320-280}{\sigma}\right)\approx0.9512. $$

但

$$ \begin{aligned}{\varPhi(1.6563)}&{{}\approx0.9512,}\\ {\frac{320-280}{\sigma}}&{{}\approx1.6563,}\\ \end{aligned} $$

$$ \sigma\approx\frac{320-280}{1.6563}=24.150. $$

我们估算一下这个考生名次位置。由于

$$ \begin{aligned}{P(X\geqslant333)}&{{}=1-P(X\lt 333)=1-\varPhi\left(\frac{333-280}{\sigma}\right)}\\ {}&{{}=1-\varPhi(2.195)\approx1-0.9859=0.0141,}\\ \end{aligned} $$

表明高于这名考生的考生所占比率大约是0.0141，在512名考生中，排名在这个考生前的人数大约为 $ 512 \times 0.0141 = 7.2167 $，即该考生大约在第9名。由此，该考生有可能被录取为免费生。

R 代码：

ex9<-function(k,m,s,x0,x1){

d<-(x0-m)/qnorm(k/s,lower.tail=FALSE)

#d<-(x0-m)/qnorm(k/s,mean=0,sd=1,lower.tail=FALSE)

rx1<-s*(1-pnorm((x1-m)/d))

rx2<-floor(rx1)+1;rx2}

ex9(25,280,512,320,333) #k<-25 m<-280 s<-512 x0<-320 x1<-333

例 9.2.10 设工厂生产的一种元件其寿命 X (单位:小时) 的密度函数为

$$ f(x)=\{\begin{aligned}&A\exp\{-\frac{x}{100}\},&&x\gt 0,\\ &0,&& 否则 .\end{aligned}. $$

求：

(1) 常数 A.

(2) 该种元件的寿命介于 50 到 150 小时之间的概率.

(3) 该种元件寿命不超过 120 小时的概率.

解 (1) $ 1 = \int_{-\infty}^{\infty} f(x) dx = \int_{0}^{\infty} A e^{-\frac{x}{100}} dx = 100A,\quad A = 0.01. $

(2) $ P(50 \lt X \lt 150) = \int_{50}^{150} \frac{1}{100} e^{-\frac{x}{100}} dx = 0.3834. $

(3) $ P(X \lt 120) = \int_{0}^{120} \frac{1}{100} e^{-\frac{x}{100}} dx = 0.00995. $

R 代码：

fx10<-function(x){

1/100*exp(-x/100)*(x>0)

}

p1<-integrate(fx10,lower=50,upper=150)

p2<-integrate(fx10,lower=0,upper=1)

list(p1=p1,p2=p2)

例 9.2.11 将一个 1 米长的木条随意截为两段，问较短的一段平均多长？

解 记 X 为截断点到左端点的距离, 可以认为 $ X \sim U(0,1) $, 本题即求

$$E(\min(X,1-X))=\int_{0}^{0.5}x\mathrm{d}x+\int_{0.5}^{1}(1-x)\mathrm{d}x=\frac{1}{8}+\frac{1}{8}=\frac{1}{4}.$$

R 代码：

fx11<-function(x){

ifelse (x<0.5,x,1-x)

y<-numeric(length(x))

y[x>0.5]<>-1-x[x>0.5]

y[x<=0.5]<>-x[x<=0.5]

y

}

curve(fx11,0,1)

p11<-integrate(fx11,lower=0,upper=1);p11 #数值解

ex11<-function(n=10000) #模拟

{

x<-runif(n)

y<-numeric(length(x))

y[x>0.5]<-1-x[x>0.5]

y[x<=0.5]<-x[x<=0.5]

sum(y)/n

}

ex11()

例 9.2.12 某日用方便食品在长时间促销期内每100个封闭的包装中装有一枚36计单个计谋卡片，宣称集齐一整套36张卡片的用户将获得大奖一份。问用户靠自己消费集齐该套卡片平均需买多少件这个公司的方便食品？

解 记 n 为整套不同卡片个数， $ X_{i} $ 为用户从集到 i-1 种不同卡片起始，到集到第 i 种卡片（不同于前面已收到的 i-1 种的任一种）时购买该种包装食品的个数，则在这期间的每次收集到的是不同于前面已收到的那 i-1 种的新卡片的概率是 $ p_{i} = \frac{n-i+1}{n} $。首次收集到不同于前面那 i-1 种的新一种卡片所用的次数即为 $ X_{i} $，服从成功率为 $ p_{i} $ 几何概率，故 $ E[X_{i}] = \frac{1}{p_{i}} $，而收集到全部 n 张卡片需购买该种包装食品的个数 $ S_{n} $ 可以表为

$$ \begin{aligned}S_{n}&=\sum_{i=1}^{n}X_{i},\\E[S_{n}]&=\sum_{i=1}^{n}E[X_{i}]=\sum_{i=1}^{n}\frac{n}{n-i+1}=n\left(1+\frac{1}{2}+\cdots+\frac{1}{n}\right).\end{aligned} $$

易知 $ E[S_n] \to \infty $。若用户仅凭自己购物得奖，对 n=36，平均大约需要购 $ 100 \times E[S_{36}] \approx 15029 $ 个该种食品。

R 代码：

ex12<-function(n=36){ # n<-36

x<-1:n

k<-n*sum(1/x)

m<-floor(100*k)+1;m}

ex12()

例 9.2.13 设 $ \{X_{i}\} $ 独立同分布, $ X_{i} \sim b(1, p) $. 将 $ X_{1}X_{2} \cdots X_{n}X_{n+1} $ 观测列出, 例如: 110100011110011 是一个可能的结果, 我们称这个顺序串中 11, 0, 1, 000, 1111, 00, 11 为游程, 这个串中有 7 个游程. 记

$$ Y_{i}=\{\begin{aligned}&1,&X_{i}\ 与 \ X_{i+1}\ 相异 ,\\ &0,& 否则 ,\end{aligned}. $$

$ S=\sum_{i=1}^{n}Y_{i} $ 为这个串 $ \{X_{i}\} $ 中出现状态变化的时点的个数，求 E[S] 和 Var[S]

解

$$ \begin{aligned}P(Y_{i}=1)&=P(X_{i}=0)P(Y_{i}=1|X_{i}=0)+P(X_{i}=1)P(Y_{i}=1|X_{i}=1)\\&=qP(X_{i+1}=1|X_{i}=0)+pP(X_{i+1}=0|X_{i}=1)\\&=qP(X_{i+1}=1)+pP(X_{i+1}=0)\\&=2pq.\\ \end{aligned} $$

$ Y_{i}\sim b(1,2pq),|j-i|\gt 1 $ 时， $ Y_{i} $ 与 $ Y_{j} $ 独立。 $ Y_{i}^{2}=Y_{i}\sim b(1,2pq) $

$$ \begin{aligned}P(Y_{i+1}Y_{i}=1)&=P(Y_{i+1}=1,Y_{i}=1)\\&=P(X_{i}=0)P(Y_{i}=1,Y_{i+1}=1|X_{i}=0)+\\&\quad P(X_{i}=1)P(Y_{i}=1,Y_{i+1}=1|X_{i}=1)\\&=P(X_{i}=0)P(X_{i+1}=1,X_{i+2}=0|X_{i}=0)+\\&\quad P(X_{i}=1)P(X_{i+1}=0,X_{i+2}=1|X_{i}=1)\\&=P(X_{i}=0)P(X_{i+1}=1)P(X_{i+2}=0)+\\&\quad P(X_{i}=1)P(X_{i+1}=0)P(X_{i+2}=1)\\&=qpq+pqp=pq.\end{aligned} $$

$$ Y_{i}Y_{i+1}\sim b(1,p q),E[S]=E\left[\sum_{i=1}^{n}Y_{i}\right]=\sum_{i=1}^{n}E[Y_{i}]=2n p q=2n p(1-p), $$

$$ \begin{aligned}E[S^{2}]&=E\left[\sum_{i=1}^{n}Y_{i}\right]^{2}=\sum_{i=1}^{n}E[Y_{i}^{2}]+2\sum_{i=1}^{n-1}E[Y_{i}Y_{i+1}]+\sum_{|j-i|\gt 1}E[Y_{i}Y_{j}]\\&=2npq+2(n-1)pq+[n^{2}-n-2(n-1)]E[Y_{i}]E[Y_{j}]\\&=2pq(2n-1)pq+4(n^{2}-3n+2)p^{2}q^{2}.\\ \end{aligned} $$

$$ \begin{aligned}\operatorname{Var}[S]&=E[S^{2}]-E^{2}[S]\\&=2pq(2n-1)pq+4(n^{2}-3n+2)p^{2}q^{2}-4n^{2}p^{2}q^{2}\\&=4npq(1-3pq)-2pq+8p^{2}q^{2}\\&=(4n-2)p-(16n-10)p^{2}+(24n-16)p^{3}-(12n-8)p^{4}.\end{aligned} $$

R 代码：

ex13<-function(m=10000,n,p){

q=1-p

s<-rep(0,m)

for (j in 1:m){

x<-sample(0:1,n+1,replace=TRUE,prob=c(p,q))

w<-rep(0,n)

w<-x[-1]-x[-(n+1)]

s[j]<sum((w!=0))

}

list(s.mean 模拟 = mean(s), s.Var 模拟 = var(s), Es=2*n*p*q, Vars= 4*n*p*q*(1-3*p*q)-2*p*q+8*p^2*q^2)

}

ex13(30000,31,1/3) # p<-1/3

例 9.2.14 为对某地 N 个人进行某种疾病调查，该种疾病可通过血样分析，分析的结果有阴性和阳性。已知根据以往经验知该地该种疾病的发病率为 p。现为加快速度拟混组验血：将受验者的血样分组检查，拟 k 个人一组，每个组将受验者的血样一半混合成该组血样，然后检查这个混合血样，如呈阳性则将该组每个人的另外一半血样逐个再检验。这种方法何时能减少工作量？

解 记 X 为每个人需验血的次数，则

$$ X\sim\left(\begin{matrix}{\displaystyle\frac{1}{k}}&{\displaystyle1+\frac{1}{k}}\\ {}&{}\\ {(1-p)^{k}}&{1-(1-p)^{k}}\\ \end{matrix}\right). $$

平均验血次数

$$ E(X)=\frac{1}{k}\times(1-p)^{k}+\left(1+\frac{1}{k}\right)[1-(1-p)^{k}]=1-(1-p)^{k}+\frac{1}{k}. $$

若 $ E(X)=1-(1-p)^{k}+\frac{1}{k}\lt 1 $ ，即 $ (1-p)^{k}\gt \frac{1}{k} $ 时，确实可以有更少的检验次数. 为此取

$$ K=\arg\min\{1-(1-p)^{k}+\frac{1}{k}\} $$

时，平均验血次数最少. 例如，对 p=0.01，取 k=11，可减少 80% 的工作量. □

R 代码：

blood<-function(p=0.1, k)

{1-(1-p)^k+1/k}

test.blood<-function(p=0.1, n=100) #n<-100

{

k<-1:n

plot(k,blood(p,k), xlim=c(1,n), ylim=c(0,2), xxt="n", xlab="组中成员数", ylab="成员平均检测次数")

axis(1,at=seq(0,n,5), labels=as.character(seq(0,n,5)))

lines(k,blood(p,k))

list(L.argmin=which.min(Blood(p,k)))

}

test.blood(0.01,100)

对 p=0.01，运行结果中的图如图 9.4 所示，表明取 k=11 时最优.

例 9.2.15 设随机向量 $ (X,Y) $ 服从 $ [-1,1]\times[-1,1] $ 上的均匀分布，求

(1) $ P(X^{2}+Y^{2}\leq1) $.

(2) 求 $ P(X^{2}+Y^{2}\leq1|Y=y) $.

(3) 验证 $ E[\sqrt{1-Y^{2}}]=P(X^{2}+Y^{2}\leqslant1) $.

解 (1) $ p(x,y)=\{\begin{aligned}&\frac{1}{4}&(x,y)\in[-1,1]\times[-1,1],\\&0& \text{其他},\end{aligned}. $ 可见 X,Y 独立.

$$ P(X^{2}+Y^{2}\leqslant1)=E(1_{X^{2}+Y^{2}\leqslant1})=\iint\limits_{X^{2}+Y^{2}\leqslant1}\frac{1}{4}\mathrm{d}x\mathrm{d}y=\frac{\pi}{4}. $$

$$ \begin{aligned}(2)P(X^{2}+Y^{2}&\leqslant1|Y=y)=P(-\sqrt{1-y^{2}}\leqslant X\leqslant\sqrt{1-y^{2}})=\int_{-\sqrt{1-y^{2}}}^{\sqrt{1-y^{2}}}\frac{1}{2}\mathrm{d}x=\\ \sqrt{1-y^{2}}.\end{aligned} $$

$$ \begin{aligned}E[\sqrt{1-Y^{2}}]&=\int_{-1}^{1}\sqrt{1-y^{2}}\times\frac{1}{2}\mathrm{d}y=\int_{0}^{1}\sqrt{1-y^{2}}\mathrm{d}y\\&=\int_{0}^{\pi/2}\cos t\cdot\cos t\mathrm{d}t=\left[\frac{t+\frac{1}{2}\sin(2t)}{2}\right]_{0}^{\pi/2}=\frac{\pi}{4}.\end{aligned} $$

R 代码：

set.seed(12345)

example.Pi1<-function(n=10000){ #模拟方法算圆周率近似值

u<-runif(n)

v<-runif(n)

x<-2*u-1

y<-2*v-1

4*sum(x^2+y^2<1)/n

}

example.Pi1(3000)

example.Pi2<-function(n=10000){ #模拟方法算圆周率近似值

v<-runif(n)

y<-2*v-1

4*sum(sqrt(1-y^2))/n

}

example.Pi2(3000)

例 9.2.16 假设一条生产线生产的产品合格率是 0.8. 要使不小于 90% 的概率保证抽验该生产线的一批产品其合格率达到在 76% 与 84% 之间，问这批产品至少要多少件？试用如下两种指定的方法求解. (1) 使用切比雪夫不等式. (2) 使用中心极限定理.

解 $ \{X_{i}\} $ 独立同分布. $ X_{i}\sim B(1,p),p=0.8,\overline{X}=\frac{1}{n}\sum_{i=1}^{n}X_{i} $

(1) 使用切比雪夫不等式.

$$ \begin{aligned}P(0.76\lt \overline{X}\lt 0.84)&=P(|\overline{X}-p|\lt 0.04)\geqslant1-\frac{\mathrm{Var}(\overline{X})}{0.04^{2}}\\&=1-\frac{p(1-p)}{0.0016n}=1-\frac{100}{n},\end{aligned} $$

为使 $ P(0.76 \lt \overline{X} \lt 0.84) \geq 0.9 $，只要 $ 1 - \frac{100}{n} \geq 0.9 $，故 $ n \geq 1000 $.

(2) 使用中心极限定理.

$$ \begin{align*}P(0.76\lt \overline{X}\lt 0.84)&=P(|\overline{X}-p|\lt 0.04)\\&=P\left(\frac{\sqrt{n}|\overline{X}-p|}{\sqrt{p(1-p)}}\leqslant\frac{0.04\sqrt{n}}{\sqrt{p(1-p)}}\right)\approx2\varPhi\left(\frac{0.04\sqrt{n}}{\sqrt{p(1-p)}}\right)-1,\end{align*} $$

为使 $ P(0.76 \lt \overline{X} \lt 0.84) \geq 0.9 $，近似地，只要 $ 2\Phi\left(\frac{0.04\sqrt{n}}{\sqrt{p(1-p)}}\right) - 1 \geq 0.9 $，即 $ \Phi\left(\frac{0.04\sqrt{n}}{\sqrt{p(1-p)}}\right) \geq 0.95 = \Phi(1.645) $，所以 $ \frac{0.04\sqrt{n}}{\sqrt{p(1-p)}} \geq 1.645 $，于是 $ n \geq \frac{1.645^{2}}{0.04^{2}} p(1-p) = \frac{1.645^{2}}{0.04^{2}} \times 0.8 \times 0.2 \approx 270.6 $，故 $ n \geq 271 $。

R 代码：

ex16.1<-function(p=0.8,e=0.04,a=0.9){

n<-p*(1-p)/(e^2*(1-a));n}

ex16.1()

ex16.2<-function(p=0.8,e=0.04,a=0.9){

s<-(1+a)/2

n<-p*(1-p)*qnorm(s,mean=0,sd=1)^2/e^2;n}

ex16.2()

例 9.2.17 某保险公司接受了 10000 辆电动自行车的保险，每辆每年的保费为 12 元。若车丢失，则车主得赔偿 1000 元。假设车的丢失率为 0.006，对于此项业务，试利用中心极限定理，求

(1) 保险公司亏损的概率 $ \alpha $.

(2) 保险公司一年获利润不少于 40000 元的概率 $ \beta $.

(3) 保险公司一年获利润不少于 60000 元的概率 $ \gamma $

解 设 X 为需要赔偿的车主人数，则需要赔偿的金额为 Y = 0.1X (万元). 保费总收入 c=12 万元. 易见，随机变量 X 服从参数为 $ (n,p) $ 的二项分布，其中 n=10000, p=0.006. $ E[X]=np=60 $, $ \mathrm{Var}(X)=np(1-p)=59.64 $. 由棣莫弗 - 拉普拉斯定理知，随机变量 X 近似服从正态分布 N(60, 59.64); 随机变量 Y 近似服从正态分布 N(6, 0.5964).

(1) 保险公司亏损的概率.

$$ \begin{aligned}\alpha&=P(Y\gt 12)=P\left(\frac{Y-6}{\sqrt{0.5964}}\gt \frac{12-6}{\sqrt{0.5964}}\right)\\&=P\left(\frac{Y-6}{\sqrt{0.5964}}\gt 7.77\right)=1-\Phi(7.77)\approx0.\end{aligned} $$

(2) 保险公司一年获利润不少于 4 万元的概率.

$$ \beta=P(12-Y\geqslant4)=P(Y\leqslant8)=P\left(\frac{Y-6}{\sqrt{0.5964}}\leqslant\frac{8-6}{\sqrt{0.5964}}\right)\approx\Phi(2.59)=0.9952. $$

(3) 保险公司一年获利润不少于 6 万元的概率.

$ \gamma = P(12 - Y \geqslant 6) = P(Y \leqslant 6) = P\left(\frac{Y - 6}{\sqrt{0.5964}} \leqslant 0\right) \approx \Phi(0) = 0.5. $

R 代码：

#int.k<function(n=10000, s=12, t=0.1, p=0.006, L1=4, L2=6){

#用中心极限定理近似计算

u<-n*p*t

v<-n*p*(1-p)*t^2

a1<1-pnorm((s-u)/sqrt(v))

a2<-pnorm((s-L1-u)/sqrt(v))

a3<-pnorm((s-L2-u)/sqrt(v))

list(A=a1, B=a2, C=a3)

}

ex17.1()

ex17.2<-function(n=10000, s=12, t=0.1, p=0.006, L1=4, L2=6){

#直接用二项分布计算

u<-p

a1<-1-pbinom((s-1)/t,n,u)

a2<-pbinom((s-L1)/t,n,u)

a3<-pbinom((s-L2)/t,n,u)

list(A=a1, B=a2, C=a3)

}

ex17.2()

例 9.2.18 (敏感性调查问题的 Simmons 模型) 涉及隐私和个人利益等方面的公众调查，直接询问受访者这类问题，受访者往往会有意答错或拒答，调查者难以收集到真实情况的数据。对这类问题的研究，最早是 Warner 提出了一种随机化问答的方法。这里描述的 Simmons 模型仍是基于 Warner 随机化问答方法的思想。

以调查学生考试作弊的比率 p 为例，设计两个问题卡片，例如

卡片 A: 你在考试中作弊了吗?

卡片 B: 你父亲的生日是在偶数月份吗?

受访者回答（限于两种结果：“是”或“否”）：[]（填写“Y”或“N”）.

其中 A 是要调查的问题，B 是为保护隐私而设计的一个与 A 无关的干扰问题，并且问题 B 回答是 “是” 的概率是知道的。要求调查者准备一个随机化装置，例如一个装有红、黑两种颜色小球的盒子，其中盒中红球的比率可以预先设定为 r。调查时给受访者展示盒中有两种同质但色标不同的球，摇匀盒中球，告

诉受访者按其抽球的颜色不同在卡片上填答一个标记“Y”或“N”，其中抽中红球仅针对问题A填答，而抽中黑球仅针对问题B填答。让受访者在相对隔离的空间抽球、并窥视其颜色，填答卡片。这使得调查者并不知道受访者回答的“Y”或“N”究竟是关于A还是关于B的。为讨论问题一般化，假定设计的回答问题B是“Y”的概率是已知的值s，根据前面设定的卡片中可近似认为 $ s\approx1/2 $。假定受访者是随机选取的可以近似作为简单随机样本，并且能在采用上述描述的方法后消除顾虑而如实填答。

(1) 求 p 的最大似然估计.

(2) 求 p 的最大似然估计 $ \hat{p}_{M} $ 的期望和方差.

(3) 利用极限定理求 p 的一个置信度为 $ 1 - \alpha $ 的近似置信区间.

(4) 如果要求 $ \hat{p}_{M} $ 的方差不超过 b，确定满足这个精度要求至少所需的样本容量.

解 若第 i 个卡片填答的是 “Y”，定义 $ X_{i}=1 $ 。否则 $ X_{i} $ 是 $ 0.X_{i}\sim b(1,\theta) $，且相互独立。 $ \theta=P(X_{i}=1) $。则由全概率公式知 $ \theta=rp+(1-r)s $。

(1) 易见 $ \theta = P(X_{i} = 1) $ 的最大似然估计是 $ \overline{X} = \frac{1}{n}(X_{1} + X_{2} + \cdots + X_{n}) $，故由最大似然估计的不变性知 $ \hat{p}_{M} = \frac{\overline{X} - (1 - r)s}{n} $ 是 p 的最大似然估计.

(2) $ E[\hat{p}_{M}] = \frac{E[\overline{X}] - (1 - r)s}{r} = \frac{\theta - (1 - r)s}{r} = p $，这表明 $ \hat{p}_{M} $ 是 p 的无偏估计。且

$$ \mathrm{Var}[\hat{p}_{M}]=\frac{\mathrm{Var}[\overline{X}]}{r^{2}}=\frac{\theta(1-\theta)}{nr^{2}}. $$

(3) 由 Var[ $ \hat{p}_{M} $] 的表示式, 构建一个统计量

$$ S^{2}=\frac{\overline{X}(1-\overline{X})}{n r^{2}}. $$

容易验证

$$ \frac{\hat{p}_{M}-p}{S}\xrightarrow{\quad L\quad}N(0,1), $$

渐近到标准正态分布，故近似地有

$$ P\left(\left|\frac{\hat{p}_{M}-p}{S}\right|\leqslant u_{\frac{\alpha}{2}}\right)\approx1-\alpha. $$

于是，可给出 p 的一个置信水平 $ 1-\alpha $ 的近似置信区间 $ (\hat{p}_{M}-u_{\frac{\alpha}{2}}S,\hat{p}_{M}+u_{\frac{\alpha}{2}}S) $

其中 $ S^{2}=\frac{\overline{X}(1-\overline{X})}{nr^{2}} $

(4) 由于 $ \mathrm{Var}(\hat{p}_{M})=\frac{\theta(1-\theta)}{nr^{2}}\leq\frac{1}{4nr^{2}} $，只要使 $ \frac{1}{4nr^{2}}\leq b $，即 $ n\geq\frac{1}{4br^{2}} $ 即可.

R 代码：

ex.Warner<-function(m,n,r,s,a=0.05){

pM<-(m/n-(1-r)*s)/r

w<-qnorm(a/2, lower.tail=FALSE)

S<-sqrt(m/n*(1-m/n)/(n*r^2))

pl<-pM-w*S

pu<-pM+w*S

list(pm=pM, pm.conf=c(pl,pu))

}

ex.Warner(90,300,0.8,0.5,0.05)

ex.n<-function(b,r){

n<-floor(1/(4*b*r^2))+1

list(n.min=n)

}

ex.n(0.05,0.5)

例 9.2.19 某人从居住地到上班地点乘坐公共交通有两种方式，第一种选择搭乘公共汽车，中途无需换乘，要穿过闹市区，路程虽短，但交通拥挤。第二种选择搭乘地铁上班，中途需要换乘，路程较长，但意外阻塞较少。假设两种方式所在路上的时间均服从正态分布。该人随机记录了12次搭乘公共汽车上班所需时间（单位：分钟）：

$$ \begin{array}{l} 46,45,88,72,40,57,50,55,58,62,46,65 \end{array} $$

该人又随机记录了13次搭乘地铁上班所需时间（单位：分钟）

$$ \begin{array}{l} 61,58,55,65,59,61,68,69,59,60,63,57,67 \end{array} $$

(1) 若有 65 分钟可用，应选取哪种交通工具上班？

(2) 若有 58 分钟可用，应选取哪种交通工具上班？

解 假设此人乘公共汽车上班所花时间为 $ X \sim N(\mu_{1}, \sigma_{1}^{2}) $，此人乘坐地铁上班所花时间为 $ Y \sim N(\mu_{2}, \sigma_{2}^{2}) $，并令 $ a_{1} = P(X \leq 65), a_{2} = P(Y \leq 65), b_{1} = P(X \leq 58), a_{2} = P(Y \leq 58) $。由所给数据得到四个参数的最大似然估计为

$$ \hat{\mu}_{1}=\bar{x}=\frac{1}{12}\sum_{i=1}^{12}x_{i}=57,\quad\hat{\sigma}_{1}=\sqrt{\frac{1}{12}\sum_{i=1}^{12}(x_{i}-\bar{x})^{2}}=12.92285, $$

$$ \hat{\mu}_{2}=\bar{y}=\frac{1}{13}\sum_{i=1}^{13}y_{i}=61.69231,\quad\hat{\sigma}_{2}=\sqrt{\frac{1}{13}\sum_{i=1}^{13}(y_{i}-\bar{y})^{2}}=4.231468, $$

因为

$$ a_{1}=P(X\leqslant65)=\Phi\left(\frac{65-\mu_{1}}{\sigma_{1}}\right), $$

所以

$$ \hat{a}_{1}=\Phi\left(\frac{65-\hat{\mu}_{1}}{\hat{\sigma}_{1}}\right)=\Phi(0.6190586)\approx0.732. $$

同理有

$$ \hat{a}_{2}=\varPhi\left(\frac{65-\hat{\mu}_{2}}{\hat{\sigma}_{2}}\right)=\varPhi(0.781689)\approx0.783, $$

$$ \hat{b}_{1}=\varPhi\left(\frac{58-\hat{\mu}_{1}}{\hat{\sigma}_{1}}\right)=\varPhi(0.07738232)\approx0.531, $$

$$ \hat{b}_{2}=\varPhi\left(\frac{58-\hat{\mu}_{2}}{\hat{\sigma}_{2}}\right)=\varPhi(-0.872583)\approx0.191. $$

从计算结果看，若有65分钟可用，应乘地铁上班。若有58分钟可用，应乘公交车上班。

R 代码：

x1<-c(46,45,88,72,40,57,50,55,58,62,46,65)

x2<-c(61,58,55,65,59,61,68,69,59,60,63,57,67)

x1.bar<-mean(x1)

x2.bar<-mean(x2)

x1.sigm<-sqrt(sum((x1-x1.bar)^2)/length(x1)) #采用方差的最大似然估计

x2.sigm<-sqrt(sum((x2-x2.bar)^2)/length(x2)) #采用方差的最大似然估计

for (xx in c(65,58))

cat("时间限制",xx," 时 P(x<="x,x",") 的估计：

乘公汽",pnorm((xx-x1.bar)/x1.sigm),";

乘地铁",pnorm((xx-x2.bar)/x2.sigm),"\\n")

例 9.2.20 X 服从 $ (0,\theta) $ 上的均匀分布, $ \theta $ 是参数. $ X_{1}, X_{2}, \cdots, X_{n} $ 是来自 X 的一个样本.

(1) 求 $ \theta $ 的矩估计和最大似然估计.

(2) 证明 $ \hat{\theta}_{1}=2\overline{X} $ 和 $ \hat{\theta}_{2}=\frac{n+1}{n}X_{(n)} $ 都是 $ \theta $ 的无偏估计.

(3) $ \theta $ 的无偏估计 $ \hat{\theta}_{1}=2\overline{X} $ 和 $ \hat{\theta}_{2}=\frac{n+1}{n}X_{(n)} $ 中，哪个更有效？

解 由设知 X 的密度函数为 $ f(x,\theta)=\{\begin{aligned}&\frac{1}{\theta},&0\lt x\lt \theta,\\ &0,& 否则.\end{aligned}. $

$$ E[X]=\int_{-\infty}^{\infty}x f(x,\theta)\mathrm{d}x=\int_{0}^{\theta}\frac{1}{\theta}x\mathrm{d}x=\frac{\theta}{2}, $$

于是 $ \theta = 2E[X] $，故得 $ \theta $ 的矩估计 $ \hat{\theta}_{MME} = 2\overline{X} $.

$ \theta $ 的似然函数

$$ L(\theta)=\prod_{i=1}^{n}f(x_{i},\theta)=\{\begin{aligned}&\frac{1}{\theta^{n}},&&0\lt x_{1},x_{2},\cdots,x_{n}\lt \theta,\\ &0,&& 否则 \end{aligned}. $$

非负函数 $ L(\theta) $ 的最大值必在 $ L(\theta) \gt 0 $ 的区域达到，即在 $ \theta \gt x_1, x_2, \cdots, x_n \gt 0 $ 上达到，而在此区域上

$$ \frac{\mathrm{d}}{\mathrm{d}\theta}L(\theta)=-\frac{n}{\theta^{n+1}}\lt 0, $$

即 $ L(\theta) $ 是 $ \theta $ 的减函数，故在这个区域 $ \theta $ 越小 $ L(\theta) $ 越大，但 $ \theta \gt x_{1}, x_{2}, \cdots, x_{n} \gt 0 $，于是 $ L(\theta) $ 在 $ \theta = \max\{x_{1}, x_{2}, \cdots, x_{n}\} = x_{(n)} $ 达到最大。故 $ \theta $ 的最大似然估计是 $ \hat{\theta}_{MLE} = X_{(n)} $。

(2) $ E[\hat{\theta}_{1}] = 2E[\overline{X}] = 2E[X] = 2 \times \frac{\theta}{2} = \theta. $ 又 $ X \sim U(0, \theta) $，对 0 < x < \theta 有

$$ \begin{aligned}F_{X_{(n)}}(x)&=P(X_{n}\leqslant x)=P(X_{1}\leqslant x,X_{2}\leqslant x,\cdots,X_{n}\leqslant x)\\&=\prod_{i=1}^{n}P(X_{i}\leqslant x)=F_{X}^{n}(x)=\frac{x^{n}}{\theta^{n}},\end{aligned} $$

可知 $ f_{X_{(n)}}(x,\theta)=\{\begin{aligned}&\frac{nx^{n-1}}{\theta^{n}},&0\lt x\lt \theta,\\&0,& 否则.\end{aligned}. $

$$ E[\hat{\theta}_{2}]=\frac{n+1}{n}E[X_{(n)}]=\frac{n+1}{n}\int_{0}^{\theta}x\,f_{X_{(n)}}(x,\theta)\mathrm{d}x=\frac{n+1}{n}\int_{0}^{\theta}\frac{n x^{n}}{\theta^{n}}\mathrm{d}x=\theta. $$

于是， $ \hat{\theta}_{1}=2\overline{X} $ 和 $ \hat{\theta}_{2}=\frac{n+1}{n}X_{(n)} $ 都是 $ \theta $ 的无偏估计.

$$ \mathrm{Var}[\hat{\theta}_{1}]=\mathrm{Var}[2\overline{X}]=4\mathrm{Var}[\overline{X}]=4\times\frac{\mathrm{Var}[X]}{n}=4\times\frac{\theta^{2}}{12n}=\frac{\theta^{2}}{3n}. $$

$$ \mathrm{V a r}[\hat{\theta}_{2}]=E[\hat{\theta}_{2}^{2}]-E^{2}[\hat{\theta}_{2}]=\left(\frac{n+1}{n}\right)^{2}\int_{0}^{\theta}\frac{n x^{n+1}}{\theta^{n}}\mathrm{d}x-\theta^{2}=\frac{\theta^{2}}{n(n+1)}. $$

注意到 $ \frac{\theta^{2}}{n(n+1)} \leqslant \frac{\theta^{2}}{3n} $，知 $ \hat{\theta}_{2} $ 比 $ \hat{\theta}_{1} $ 更有效。一个模拟计算比较的示意图见图9.5. □

R 代码：

num.comp<-function(theta,k=10000,n=10){ # 模拟对方差大小的比较

y<-matrix(0,k,2)

for (i in 1:k){

x<-runif(n,0,theta)

y[i].[-c(2*mean(x),(1+n)*max(x)/n)]

c(sum((y[,1]-theta)^2),sum((y[,2]-theta)^2))/k)

对 theta 介于 1~10 的一些值，模拟计算 theta1 和 theta2 的方差，作图显示

set.seed(123)

theta<-seq(1,10,by=0.2); s<-length(theta); yy<-matrix(0,s,2)

for (j in 1:s){

yy[j,1]<-num.comp(theta[j],) [1]

yy[j,2]<-num.comp(theta[j],) [2]

plot(theta, yy[,1],type="b",ylim=c(0,3),pch=1,1ty=1,

xlab= expression(paste(theta,"的值")),ylab="估计量的方差")

lines(theta, yy[,2],type="b",1ty=1,pch=2)

ex1<-expression(paste("Var",hat(theta)[1],")")

ex2<-expression(paste("Var",hat(theta)[2],")")

legend(x=1,y=2.9,legend=c(ex1,ex2),pch=c(1,2))

例 9.2.21 设总体服从尺度参数为 1 的柯西分布，即总体的密度函数为 $ p(x)=\frac{1}{\pi[1+(x-p)^{2}]} $，试模拟求位置参数 p 的最大似然估计.

解 由于柯西总体位置参数 p 的最大似然估计没有解析表达，下面以柯西总体位置参数 p = 1 的情形为例进行抽样，并对所抽样本用最大似然估计计算估计值。我们两种方式求最大似然估计值，分别是：(1) 求似然方程根的数值解。

(2) 求似然函数的最大值点的数值解.

(1) 求似然方程根的数值解.

利用下面的 R 程序:

利用下面的 R 程序:

set.seed(123)

x<-rcauchy(100,1)

fx20<-function(p) sum((x-p)/(1+(x-p)^2))

out<-uniroot(f,c(0,5)) #求 f 在 (0,5) 上的根

out

结果显示:

root

[1] 1.040760

f.root

[1] 7.958795e-07

iter

[1] 6

estim.prec

[1] 6.103516e-05

这次取样本容量 n = 100, 得到最大似然估计值 $ \hat{p}_1 = 1.040760 $.

(2) 求似然函数的最大值点的数值解

利用下面的 R 程序:

set.seed(123)

x<-rcauchy(100,1)

loglike<-function(p)-sum(log(1+(x-p)^2)) #构建对数似然函数，略去了常数项

optimize(loglike,c(0,5), maximum = TRUE) #求对数似然函数在 (0,5) 上最大点

结果显示:

maximum

[1] 1.040754

objective

[1] -129.7924

得到最大似然估计值 $ \hat{p}_{2}=1.040754 $

注：为了方便对照，我们指定了同一种子的随机数。两种不同方法得到的 p 的估计值与选定的 p 的真值 1 很接近。

例 9.2.22 设总体 $ X \sim N(\mu,1) $. 试模拟用来自于 $ N(\mu,1) $ 的正态总体样本, 计算 $ \mu $ 的度为 0.95 的置信区间, 以理解置信区间与置信度的意义.

解 假设 $ X_{1}, X_{2}, \cdots, X_{n} $ 是来自总体 X 的样本，则 $ \mu $ 的置信度为 0.95 的置信区间为

$$ \left[\overline{{X}}-\frac{u_{0.975}}{\sqrt{n}},\overline{{X}}+\frac{u_{0.975}}{\sqrt{n}}\right]=\left[\overline{{X}}-\frac{1.96}{\sqrt{n}},\overline{{X}}+\frac{u_{0.975}}{\sqrt{n}}\right], $$

即 $ P\left(\overline{X}-\frac{u_{0.975}}{\sqrt{n}}\lt \mu\lt \overline{X}+\frac{u_{0.975}}{\sqrt{n}}\right)=0.95. $

该置信区间会随着样本观测值的不同而不同，但能包含 $ \mu $ 的置信区间的概率为 0.95.

为方便起见，我们从总体 $ N(0,1) $ 中（即 $ \mu = 0 $）随机生成 200 个样本观测值，由此估计总体的均值和均值的 95% 置信区间。把该试验重复实验 100 次。请阅读本题所附的程序代码。为便于读者对照，我们指定了种子。从运行结果看，这 100 次实验所得估计值与真值误差平方的平均值是 0.003220607。而得到的 100 个置信区间有 94 次包含了 $ \mu $ 的均值 0，只有 6 次没有包含 $ \mu $ 的真值。在图 9.6 中，横坐标表示试验的序数，纵线表示置信区间线段。

R 代码：

interval.graph<-function(xl,xu){

y1=xl;y2=xu;n=length(y1)

plot(y1,type="n",ylim=c(-.3,.3),xlab="",ylab="")

segments((1:n)[y1<0&y2>0],y1[y1<0&y2>0],(1:n)[y1<0&y2>0],y2[y1<0&y2>0],col="blue")

segments((1:n)[y1>0],y1[y1>0],(1:n)[y1>0],y2[y1>0],col="red",lwd=3)

segments((1:n)[y2<0],y1[y2<0],(1:n)[y2<0],y2[y2<0],col="red",lwd=3)

SUM=sum(xl<=0 & xu>=0)

abline(h=0)

cat("这",n,"次实验所得95%区间估计包含真值的个数是",SUM,"\n")}

set.seed(333)

n=200 #每次实验的样本容量

m=100 #实验次数

x=rep(0,m) #x=vector("numeric",10) 或 x=array(0,m)

xl=rep(0,m)

xu=rep(0,m)

for (i in 1:m){

x[i]=mean(rnorm(n))

xl[i]=x[i]+qnorm(0.025)*sqrt(1/n)

xu[i]=x[i]+qnorm(0.975)*sqrt(1/n)

}

err.squ=sum(x^2)/n

cat("这",n,"次实验所得估计值与真值误差平方的平均是",err.squ,"\n")

interval.graph(xl,xu)

例 9.2.23 从某小学五年级男生中抽取 72 人，测量身高，得数据 (单位：cm) 如下：

128.1, 144.4, 150.3, 146.2, 140.6, 139.7, 134.1, 124.3, 147.9, 143.0, 143.1, 142.7, 126.0, 125.6, 127.7, 154.4, 142.7, 141.2, 133.4, 131.0, 125.4, 130.3, 146.3, 146.8, 142.7, 137.6, 136.9, 122.7, 131.8, 147.7, 135.8, 134.8, 139.1, 139.0, 132.3, 134.7, 150.4, 142.7, 144.3, 136.4, 134.5, 157.3, 152.7, 148.1, 139.6, 138.9, 136.1, 135.9, 142.2, 152.1, 142.4, 142.7, 136.2, 135.0, 154.3, 147.9, 141.3, 143.8, 138.1, 139.7, 127.4, 146.0, 155.8, 141.2, 146.4, 139.4, 140.8, 127.7, 150.7, 160.3, 148.5, 162.5.

按以往经验，五年级男生身高近似服从正态分布，现要从五年级男生中挑选身高在143至153的组成健美操方队。

(1) 试估计在五年级男生中任选一人身高在 145 至 153 间的概率.

(2) 求五年级男生平均身高 $ \mu $ 的置信度为 0.95 的置信区间.

(3) 五年级男生身高是否服从正态分布？

解 (1) 设任从五年级男生挑选一男生身高为 X，并设 $ X \sim N(\mu, \sigma^{2}) $，记

$$ p=P\{145\leqslant X\leqslant153\}, $$

则 $ p = \Phi\left(\frac{153 - \mu}{\sigma}\right) - \Phi\left(\frac{145 - \mu}{\sigma}\right) $，又

$$ \bar{x}=\frac{1}{72}\sum_{i=1}^{72}x_{i}=140.6889,\quad s=\sqrt{\frac{1}{72}\sum_{i=1}^{72}(x_{i}-\bar{x})^{2}}=8.7270, $$

则 $ \mu,\sigma $ 的最大似然估计值分别为

$$ \hat{\mu}=\bar{x}=140.6889,\quad\hat{\sigma}=s=8.7270, $$

$$ \hat{p}=p=\Phi\left(\frac{153-\mu}{\sigma}\right)-\Phi\left(\frac{145-\mu}{\sigma}\right)=\Phi(1.411)-\Phi(0.494)=0.234. $$

(2) $ \mu $的置信度为0.95的置信区间为

$$ \left[\bar{x}-\frac{s^{*}}{\sqrt{n}}t_{0.975}(71),\bar{x}+\frac{s^{*}}{\sqrt{n}}t_{0.975}(71)\right]=[138,142]. $$

(3) 用 $ F_{0}(x) $ 表示正态分布 $ N(\mu,\sigma^{2}) $ 的分布函数，则本题检验的假设问题为 $ H_{0} $ ：五年级男生身高的分布函数为 $ F_{0}(x) $， $ H_{1} $ ：五年级男生身高的分布函数不是 $ F_{0}(x) $ 。而参数 $ \mu,\sigma $ 的最大似然估计为

$$ \hat{\mu}=\bar{x}=140.6889,\quad\hat{\sigma}=s=8.7270. $$

在 $ N(140.6889,8.7270^{2}) $ 分布下，给出总体落在区间 $ (a_{i-1},a_{i}] $ 内概率的估计值

$$ \hat{p}_{i}=\varPhi\left(\frac{a_{i}-140.6889}{8.7270}\right)-\varPhi\left(\frac{a_{i-1}-140.6889}{8.7270}\right). $$

计算结果如表9.3所示.

从表 9.3 可知 $ \chi^{2} $ 检验统计量观测值为 0.65629, p 值 = 0.88 > 0.05, 即不拒绝正态性假设. □

R 代码：

data.shg<-

c (128.1,144.4,150.3,146.2,140.6,139.7,134.1,124.3,147.9,143.0,

143.1,142.7,126.0,125.6,127.7,154.4,142.7,141.2,133.4,131.0,

125.4,130.3,146.3,146.8,142.7,137.6,136.9,122.7,131.8,147.7,

135.8,134.8,139.1,139.0,132.3,134.7,150.4,142.7,144.3,136.4,

134.5,157.3,152.7,148.1,139.6,138.9,136.1,135.9,142.2,152.1,

142.4,142.7,136.2,135.0,154.3,147.9,141.3,143.8,138.1,139.7,

127.4,146.0,155.8,141.2,146.4,139.4,140.8,127.7,150.7,160.3,

148.5,162.5)

n<-length(data.shg)

alpha<-0.05

x.mean<-mean(data.shg); x.sd<-sd(data.shg)

xx<-c(145,153)

cat("身高介于 145-153 间的概率估计：",

pnorm((xx[2]-x.mean)/x.sd)-pnorm((xx[1]-x.mean)/x.sd),"\n")

T<-qt(1-alpha/2,n-1)

rr<-round(c(x.mean-x.sd*T/sqrt(n),x.mean+x.sd*T/sqrt(n)),2)

data.frame(置信下限 =rr[1], 置信上限 =rr[2], row.names='置信区间')

#或如下直接用 t.test() 给出置信区间

t.test(data.shg,conf.level=1-alpha)$conf$

r<-2 #F_0 中参数个数

b<-c(-Inf,130,135,140,145,150,Inf)

a<-table(cut(data.shg,br=b))

p<-pnorm(b[-1],mean(data.shg),sd(data.shg)*sqrt((n-1)/n))

p.hat<-p-c(0,p[-length(p)])

m<-length(a)

a.chi<-sum((a-n*p.hat)^2/(n*p.hat))

a.p<-pchisq(a.chi,df=m-r-1, lower.tail=FALSE)

data.frame(x^2检验统计值=a.chi,样容量=n,自由度=m-r-1,

p值=a.p,row.names="%)

结果显示：

身高介于 145-153 间的概率估计：0.2312430

置信下限 置信上限

置信区间 138.62 142.75

$$ x^{2} $$

p 值=0.883741，故在 $ \alpha=0.05 $ 下不拒绝正态性假设.

例 9.2.24 设总体 $ X \sim N(\mu, 2) $，对于检验问题 $ H_{0} $: $ \mu = 3 $, $ H_{1} $: $ \mu = 4 $，假定检验由拒绝域 $ C = \{(x_{1}, x_{2}, \cdots, x_{n}) : \frac{1}{n}(x_{1} + x_{2} + \cdots + x_{n}) \gt 3.8\} $ 确定.

(1) 若 n = 32，分别求检验犯第一类错误及第二类错误的概率 $ \alpha $ 和 $ \beta $.

(2) 若要求检验犯第二类错误的概率 $ \beta \leqslant 0.05 $，求 n 至少要取多少.

(3) 证明： $ n \to \infty $ 时， $ \alpha \to 0 $ 且 $ \beta \to 0 $.

解 $ H_{0}:\mu=3 $, $ H_{1}:\mu=4 $. $ C=\{(x_{1},x_{2},\cdots,x_{n}):\frac{1}{n}(x_{1}+x_{2}+\cdots+x_{n})\gt 3.8\} $.

$$ \begin{aligned}\alpha&=P( 第一类错误 )=P(\overline{X}\in C|H_{0})\\&=P\left(\frac{\overline{X}-3}{\sqrt{2}/\sqrt{n}}\gt \frac{3.8-3}{\sqrt{2}/\sqrt{n}}|X\sim N(3,2)\right)\\&=P\left(\frac{\overline{X}-3}{\sqrt{2}/\sqrt{n}}\gt 0.8\sqrt{\frac{n}{2}}|X\sim N(3,2)\right)\\&=1-\varPhi\left(0.8\sqrt{\frac{n}{2}}\right),\end{aligned} $$

$$ \begin{aligned}\beta&=P( 第二类错误 )=P(\overline{X}\notin C|H_{1})\\&=P\left(\frac{\overline{X}-4}{\sqrt{2}/\sqrt{n}}\leqslant\frac{3.8-4}{\sqrt{2}/\sqrt{n}}|X\sim N(4,2)\right)\\&=P\left(\frac{\overline{X}-4}{\sqrt{2}/\sqrt{n}}\leqslant-0.2\sqrt{\frac{n}{2}}|X\sim N(4,2)\right)\\&=\varPhi\left(-0.2\sqrt{\frac{n}{2}}\right)=1-\varPhi\left(0.2\sqrt{\frac{n}{2}}\right).\end{aligned} $$

(1) n=32 时，

$$ \alpha=1-\varPhi(3.2)\approx0.000687,\qquad\beta=1-\varPhi(0.8)\approx0.211855. $$

(2) 由设知 $ \beta = 1 - \Phi \left(0.2\sqrt{\frac{n}{2}}\right) \leqslant 0.01 $，故

$$ 0.2\sqrt{\frac{n}{2}}\geqslant\varPhi^{-1}(0.01)\approx2.326348, $$

得 $ n \geqslant 270.59 $，于是 $ n \geqslant 271 $.

(3) 由于 $ \lim_{x\to+\infty}\Phi(x)=1 $，而

$$ \alpha=1-\Phi\left(0.8\sqrt{\frac{n}{2}}\right),\quad\beta=\Phi\left(-0.2\sqrt{\frac{n}{2}}\right)=1-\Phi\left(0.2\sqrt{\frac{n}{2}}\right). $$

故证.

R 代码：

alph.beta<-function(c,n,mu1,mu2,sig.sq,beta0=0.01){

#alph1<-1-pnorm((c-mu1)*sqrt(n)/sqrt(sig.sq))

alph<-1-pnorm(c,mu1,sqrt(sig.sq/n))

beta<-pnorm(c,mu2,sqrt(sig.sq/n))

n0<-sig.sq*(qnorm(beta0))^2/(c-mu2)^2

list(第一类错误概率=alph, 第二类错误概率=beta, n0=n0, n=floor(n0)+1)}

例 9.2.25 假设某工厂有两条流水线生产同一种产品，据以往观测数据可以假定这两条流水线日产量均近似地服从同方差的正态分布，且两条流水线的运行互不干扰，现对两条流水线随机地记录了 12 天的产量（单位：件）如下：

现在要求通过这些统计数据分别求两条流水线日平均产量的置信度为95%的置信区间及两条流水线日平均产量之和的置信度为95%的置信区间。

解 假设第一条流水线日产量 $ X \sim N(\mu_{1}, \sigma^{2}) $，第二条流水线日产量 $ Y \sim N(\mu_{2}, \sigma^{2}) $，经计算得

$$ \bar{x}=498.33333,\ s_{1}^{*2}=81.51518,\ \bar{y}=503.91667,\ s_{2}^{*2}=116.26527,\ s_{w}=9.94436. $$

又 $ t_{0.975}(11)=2.2010 $, $ t_{0.975}(22)=2.5083 $, 故 $ \mu_{1} $ 的 0.95 的置信区间为

$$ \left[\bar{x}-\frac{s_{1}^{*}}{\sqrt{n}}t_{0.975}(11),\bar{x}+\frac{s_{1}^{*}}{\sqrt{n}}t_{0.975}(11)\right]=[492.596\;9,504.069\;8]. $$

$ \mu_{2} $ 的 0.95 的置信区间为

$$ \left[\bar{y}-\frac{s_{2}^{*}}{\sqrt{n}}t_{0.975}(11),\bar{y}+\frac{s_{2}^{*}}{\sqrt{n}}t_{0.975}(11)\right]=[497.0657,510.7677]. $$

$ \mu_{1}+\mu_{2} $ 的 0.95 的置信区间为

$$ \begin{aligned}{}&{{}\left[\bar{x}+\bar{y}-s_{w}t_{0.975}(22)\sqrt{\frac{1}{12}+\frac{1}{12}},\bar{x}+\bar{y}+s_{w}t_{0.975}(22)\sqrt{\frac{1}{12}+\frac{1}{12}}\right]}\\ {=}&{{}~[993.851\;7,1\;010.648\;0].}\\ \end{aligned} $$

R 代码：

x<-c(500,510,498,501,495,478,495,489,512,504,501,497)

y<-c(508,510,519,506,504,490,498,480,512,515,503,502)

t.test(x)

t.test(y) $conf.int #此处选择了只显示置信区间有关信息$

t.test(x,-y, conf.level = 0.95)$conf.int #此处选择了只显示置信区间的信息$

#或采用以下按公式自编程序求 x+y 的置信区间

conf.interv1<-function(x,y, conf.level=0.95){

n1<-length(x);n2<-length(y)

x.bar<-mean(x);y.bar<-mean(y);

xs<-var(x);ys<-var(y)

alpha<-(1+conf.level)/2

sw<-sqrt(((n1-1)*xs+(n2-1)*ys)/(n1+n2-2));sw

dd<-sw*sqrt((1/n1+1/n2))

tq<qt(alpha,n1+n2-1)

ss<-x.bar+y.bar

cat("mu1+mu2 的置信区间: ("ss-tq*dd,","ss+tq*dd,")\n")

}

x<-c(500,510,498,501,495,478,495,489,512,504,501,497)

y<-c(508,510,519,506,504,490,498,480,512,515,503,502)

conf.interv1(x,y,)

例 9.2.26 对一批由同一种纱线织成的袜子，在水温分别为 $ 30^{\circ} $C 与 $ 40^{\circ} $C 水中进行洗搓收缩率试验，其他条件完全相同，测得两种温度下袜子收缩率为

水温 $ 30^{\circ} $C: 4.3, 3.4, 6.4, 7.6, 6.5, 3.5, 6.2

水温 $ 40^{\circ} $C: 6.1, 7.3, 4.2, 4.1, 5.5

假定两种温度下袜子收缩率分别近似地服从正态分布 $ N(\mu_{1},\sigma_{1}^{2}) $， $ N(\mu_{2},\sigma_{2}^{2}) $。

(1) 试求总体方差比 $ \sigma_{1}^{2}/\sigma_{2}^{2} $ 的置信度为 95% 的置信区间.

(2) 假定两总体方差相等时，试求总体均值差 $ \mu_{1}-\mu_{2} $ 的置信度 95% 的置信区间.

解 经计算两种湿度下收缩率的样本均值与修正样本方差分别为

$$ \bar{x}=\frac{1}{7}\sum_{i=1}^{7}x_{i}=5.414286,\quad s_{1}^{*2}=\frac{1}{6}\sum_{i=1}^{7}(x_{i}-\bar{x})^{2}=2.751429, $$

$$ \bar{y}=\frac{1}{5}\sum_{i=1}^{5}y_{i}=5.440000,\quad s_{2}^{*2}=\frac{1}{4}\sum_{i=1}^{5}(y_{i}-\bar{y})^{2}=1.808000, $$

$$ s_{w}=\sqrt{\frac{6s_{1}^{*2}+4s_{2}^{*2}}{10}}=1.540798. $$

$$ F_{0.975}(6,4)=9.197,\quad F_{0.975}(4,6)=6.227,\quad t_{0.975}(10)=2.228\:1. $$

(1) $ \sigma_{1}^{2}/\sigma_{2}^{2} $ 的置信度为95%的置信区间

$$ \left[\frac{s_{1}^{*2}}{F_{0.975}(6,4)s_{2}^{*2}},\frac{s_{1}^{*2}}{s_{2}^{*2}}F_{0.975}(4,6)\right]=[0.1655\mathrm{~,~}9.4763]\mathrm{.} $$

(2) $ \mu_{1}-\mu_{2} $的置信度为95%的置信区间

$$ \begin{aligned}{}&{{}\left[\bar{x}-\bar{y}-s_{w}t_{0.975}(10)\sqrt{\frac{1}{7}+\frac{1}{5}},\bar{x}-\bar{y}+s_{w}t_{0.975}(10)\sqrt{\frac{1}{7}+\frac{1}{5}}\right]}\\ {=}&{{}\left[-2.035\;938,1.984\;510\right].}\\ \end{aligned} $$

R 代码：

$$ \begin{aligned}&x\lt -c(4.3,3.4,6.4,7.6,6.5,3.5,6.2)\\&y\lt -c(6.1,7.3,4.2,4.1,5.5)\\&var.test(x,y)\conf\quad\# 求方差比的置信区间 \\&t.test(x,y,var.equal=TRUE)\conf\quad\# 方差相等时求均值差的置信区间 \end{aligned} $$

例 9.2.27 点心厂质检员为判断牛奶供应商所提供的牛奶是否被兑水，对供应商供应的牛奶进行了随机抽样检查，测得11个鲜奶样品的冰点数据如下：

$$ \begin{aligned}&-0.543,-0.547,-0.536,-0.548,-0.544,-0.547,\\ &-0.542,-0.547,-0.549,-0.541,-0.548.\\ \end{aligned} $$

假定供应商供应的牛奶的冰点服从正态分布，并已知天然鲜奶的冰点是摄氏 -0.545 度。给定显著性水平 $ \alpha = 0.01 $。

(1) 供应商提供的牛奶是否兑水？

(2) 若以往供应商提供的牛奶其冰点标准差 $ \sigma = 0.004 $，可否认为这次所查的牛奶仍正态分布 $ N(\mu, 0.004^{2}) $?

(3) 若供应商的牛奶其冰点服从 $ N(\mu, 0.004^{2}) $，由数据判断提供的牛奶是否兑水.

解 设 $ X_{i} $ 表示第 i 个样品的冰点，可认为 $ X_{1}, X_{2}, \cdots, X_{n} $ 是来自正态总体 $ X \sim N(\mu, \sigma^{2}) $ 的样本，i = 11. 由于水的冰点是摄氏 0 度，兑水后的牛奶的冰点会有所提高.

(1) $ H_{01} $: $ \mu \leqslant \mu_{0} $, $ H_{11} $: $ \mu \gt \mu_{0} $. 这里 $ \mu_{0} = -0.545 $, $ \sigma^{2} $ 未知.

$$ t=\frac{\bar{x}-\mu_{0}}{s}=1.5301\lt t_{0.995}(10)=3.169, $$

或从 p 值 = 0.0785 > $ \alpha = 0.01 $，故不拒绝 $ H_{01} $，不认为牛奶被兑了水。

(2) $ H_{02} $: $ \sigma^{2}=\sigma_{0}^{2} $, $ H_{12} $: $ \sigma^{2}\neq\sigma_{0}^{2} $. 这里 $ \sigma_{0}^{2}=0.004^{2} $, $ \mu $ 未知.

$$ \chi^{2}=\frac{(n-1)s^{2}}{\sigma_{0}^{2}}=8.761364,\ \chi_{0.995}^{2}(10)=25.188,\ \chi_{0.005}^{2}(10)=2.156, $$

结合拒绝域的表达形式，或由 p 值 = 0.890218 > $ \alpha = 0.01 $，故不拒绝 $ H_{02} $，不认为方差有改变，认为 $ X \sim N(\mu, 0.004^{2}) $。

(3) $ H_{03} $: $ \mu \leqslant \mu_{0} $, $ H_{13} $: $ \mu \gt \mu_{0} $. 这里 $ \mu_{0} = -0.545 $, $ \sigma_{2} = 0.004^{2} $ 已知.

$$ u=\frac{\bar{x}-\mu_{0}}{\sigma/\sqrt{n}}=1.432179\lt u_{0.995}=2.576, $$

或由 p 值 = 0.07604632 > $ \alpha = 0.01 $，故不拒绝 $ H_{03} $，不认为牛奶被兑了水。

R 代码：

x<-c(-0.543,-0.547,-0.536,-0.548,-0.544,-0.547,-0.542,-0.547,-0.541,-0.541,-0.540)

alpha<-0.01

t.test(x, ,alternative="greater", mu=-0.545, conf.level=1-alpha)

n<-length(x)

sig<-0.004

s2<-var(x);df<-n-1

chisq.x<-df*s2/sig^2

p.value2<-2*min(pchisq(chisq.x, df), pchisq(chisq.x, df, lower=F))

data.frame(chisq=chisq.x, df=df, p_value=p.value2, S2=s2,

ci.L=df*s2/qchisq(1-alpha/2, df),

ci.U=df*s2/qchisq(alpha/2, df), row.names="

mu0<-0.545

mean<-mean(x)

u<-(mean-mu0)/(sig/sqrt(n))

p.value3<-1-pnorm(u)

b<-qnorm(1-alpha,mean=0, sd=1, lower=T)* sig/sqrt(n)

data.frame(u=u, p_value=p.value3, mean=mean,

ci.L=mean-b, ci.U=Inf, row.names="

例 9.2.28 火灾报警探测器的生产厂家关心的问题之一为探测器在湿度小的环境下标定的结果与在湿度大的环境下标定的结果有无显著差异. 对于这样一个在生产中急需解决的问题, 从数理统计的角度来看可归结为检验两组观测

数据的均值是否相等的问题。某厂家大批量生产离子感烟式火警器，由于在生产过程中各技术指标得到严格控制，因此探测器产品的电离电流值从理论上讲呈正态分布。表9.4和表9.5所列数据分别为火警探测器经过冲潮处理与未经冲潮处理复标结果，表中略去了数据的物理含义和计量单位。

问湿度对火灾报警探测器复标结果是否有显著影响（显著水平 $ \alpha = 0.05 $）?

解 设冲潮后的复标结果与未冲潮的复标结果分别为 X,Y, 由题目假设知 $ X \sim N(\mu_{1},\sigma_{1}^{2}) $, $ Y \sim N(\mu_{2},\sigma_{2}^{2}) $. 首先检验假设 $ H_{01}:\sigma_{1}^{2}=\sigma_{2}^{2} $, $ H_{11}:\sigma_{1}^{2} \neq \sigma_{2}^{2} $. 经计算得

$$ \bar{x}=\frac{1}{83}\sum_{i=1}^{83}x_{i}=68.3494,\quad s_{1}^{*2}=\frac{1}{82}\sum_{i=1}^{83}(x_{i}-\bar{x})^{2}=15.2301, $$

$$ \bar{y}=\frac{1}{105}\sum_{i=1}^{105}y_{i}=67.0952,\quad s_{2}^{*2}=\frac{1}{104}\sum_{i=1}^{105}(y_{i}-\bar{y})^{2}=15.5678, $$

$$ s_{w}=\sqrt{\frac{82s_{1}^{*2}+104s_{2}^{*2}}{186}}=7.3004, $$

$$ F=\frac{s_{1}^{*2}}{s_{2}^{*2}}=0.9783. $$

由于 F = 0.9783, p 值 = 0.9232 > 0.05, 故不拒绝 $ H_{01} $.

其次检验假设 $ H_{02} $: $ \mu_{1} = \mu_{2} $, $ H_{12} $: $ \mu_{1} \neq \mu_{2} $. 经计算检验统计量的值为

$$ t=\frac{\bar{x}-\bar{y}}{s_{w}\sqrt{\frac{1}{n_{1}}+\frac{1}{n_{2}}}}=2.1774, $$

p 值 = 0.031 < 0.05，故拒绝原假设 $ H_{02} $，即认为湿度对火灾报警探测器复标结果有显著影响. □

R 代码:

x1<-c(60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,78,80)

f1<-c(1,2,2,6,3,3,13,5,8,5,8,11,8,2,2,1,1,1,1)

x2<-c(58,59,60,61,62,63,64,65,66,67,68,69,70,71,73,74,75,78)

f2<-c(2,1,2,3,6,6,8,6,10,8,16,12,8,6,4,3,3,1)

x11<-rep(x1,f1)

x22<-rep(x2,f2)

var.test(x11,x22) #检验冲潮前后方差是否相同

t.test(x11,x22, var.equal = TRUE) #在方差相同假设下，检验复标结果是否相同

t.test(x11,x22) #或直接选用不做方差相等假设的 t 检验

例 9.2.29 厂家 A, B, C 是生产某一产品的知名企业，在过去的一年里，三大工厂生产的该产品的市场占有率分别为 15%，35%，25%。厂家 A 为了提高市场占有率，对该产品进行了改进，现在在市场销售厂家 A 的产品均为新型产品。现进行抽样调查，对销售出的 200 件调查结果如表 9.6。

试依据调查数据对该产品的市场占有率是否发生变化做出判断，以便为厂家 A 下一步的决策提供依据 (显著水平 $ \alpha = 0.05 $).

解 设 $ p_{1}, p_{2}, p_{3} $ 分别为厂家 A, B, C 是生产该产品在市场的占有率，在显著水平 $ \alpha = 0.05 $ 下，检验假设 $ H_{0}: p_{1} = 0.15, p_{2} = 0.35, p_{3} = 0.25, H_{1} $ ：该产品的市场占有率发生了变化。令 $ p_{4} = 1 - p_{1} - p_{2} - p_{3} = 0.25 $，依题意知 n = 210， $ v_{1} = 42, v_{2} = 67, v_{3} = 49, v_{4} = 42 $。检验统计量 $ \chi^{2} $ 的值为

$$ \chi^{2}=\sum_{i=1}^{4}\frac{(v_{i}-np_{i})^{2}}{np_{i}}=6.2286, $$

p 值 = 0.101 > 0.05，故不拒绝 $ H_{0} $，即现有数据不拒绝市场占有率未变这一论断.

R 代码：

x<-c(42,67,49,42)

pd<-c(0.15,0.35,0.25,0.25)

chisq.test(x,p=pd)

例 9.2.30 为了检验某暑期英语辅导班的效果，从某校学生中随机选取 25 名学生参加该暑期辅导班，在辅导班开始和结束时分别进行了一次难度相当的综合性考试，各学生的考试成绩见表 9.7，据以往经验知，学生参加辅导班后英语成绩之差近似服从正态分布，在显著水平 $ \alpha = 0.05 $ 下，问以下成绩能否认为参加该辅导班对英语成绩的提高是有效的？

解 设 X 表示从某校学生中随机选取一人其参加辅导班前的英语成绩, Y 表示其参加辅导班后的英语成绩, 令 Z = X - Y, $ z_{i} = x_{i} - y_{i} $ ( $ i = 1, 2, \cdots, 25 $) 为来自 Z 的一组样本观测值, 并将 $ z_{i} $ 列如下:

$$ \begin{array}{l}-2,2,-8,-7,3,-2,-8,2,-13,-1,-18,3,-1,\\ \\ -2,-3,1,-10,2,-2,2,2,-4,-3,-2,-4.\end{array} $$

现设 $ Z \sim N(\mu, \sigma^{2}) $，要检验的假设为 $ H_{0}: \mu = 0, H_{1}: \mu \lt 0 $。经计算得

$$ \bar{z}=\frac{1}{25}\sum_{i=1}^{25}z_{i}=-2.92,\quad s^{*}=\frac{1}{24}\sum_{i=1}^{25}(z_{i}-\bar{z})^{2}=5.2751. $$

检验统计量的值为 $ t = \frac{\sqrt{24} \bar{z}}{s} = -2.7677 $，p 值 = 0.00535 < $ \alpha = 0.05 $，故拒绝 $ H_{0} $，即认为参加该辅导班对英语成绩的提高有效。

R 代码：

#此处假定题目中表的数据存储在 d:\R-ex 下，是 csv 格式文件，文件名为 xuesh.csv x<-x[,-1]

$$ \operatorname{c o l n a m e s}(\mathbf{x})=\mathbf{c}(^{\prime}\mathbf{x}1^{\prime},^{\prime}\mathbf{x}2^{\prime}) $$

attach(x)

t.test(x1,x2,paired=TRUE,alternative='less')