R软件在概率论与数理统计中的应用解读word文档良心出品.docx

R软件在概率论与数理统计中的应用解读word文档良心出品.docx

《R软件在概率论与数理统计中的应用解读word文档良心出品.docx》由会员分享，可在线阅读，更多相关《R软件在概率论与数理统计中的应用解读word文档良心出品.docx（42页珍藏版）》请在冰豆网上搜索。

R软件在概率论与数理统计中的应用解读word文档良心出品

第十章R软件在概率论与数理统计中的应用

概率论与数理统计是一门应用性很强的数学学科.在实际应用中，经常会遇到比较繁琐的计算，在这种情况下，可以借助计算机编程计算.R是由奥克兰大学RossIhaka和RobertGentleman建立的统计分析软件，是AT&T贝尔实验室创建的S语言的民用版本，S-Plus是Insightful公司的商业化软件.而R是免费的，且源代码完全开放，其源代码主要用C语言编写，也有一些程序用Fortran编写，熟悉C语言的读者很容易掌握R.本章将结合概率统计教学，具体介绍R软件的使用方法.

§10.1R软件概述

10.1.1R软件的下载与安装

R软件是完全免费的，在http:

//cran.r-project.org/bin/windows/base/上可以下载到R软件的Windows版本，当前的版本是R-2.8.1版（2008年12月22日发布），大约31MB.

R软件的安装非常容易，运行刚才下载的程序，按照Windows的提示安装即可.当开始安装后，选择安装提示的语言，接受安装协议，选择安装目录，并选择安装组件.在安装组件中，最好将PDFReferenceManual项也选上，这样在R软件的帮助文件中有较为详细的PDF格式的软件说明.

安装完成后，程序会创建R程序组并在桌面上创建R主程序的快捷方式.通过快捷方式运行R软件，便可调出R软件的主窗口，如图10-1所示.

图10-1

R软件的界面与Windows的其他编程软件类似，由一些菜单和快捷按钮组成.快捷按钮下面的窗口便是命令输入窗口，它也是运算结果的输出窗口.

主窗口上方的一些文字是运行R软件时出现的一些说明和指引.文字下的“>”符号便是R软件的命令提示符，在其后可以输入命令.R软件一般采用交互式工作方式，在命令提示符后输入命令，回车后便会输出计算结果.当然也可将所有的命令建立成一个文件，运行这个文件的全部或部分来执行相应的命令，从而得到相应的结果.

对于一些陌生的命令，可以通过单击主窗口中的【帮助】菜单查找其具体的命令格式及含义.在命令的后面加一个“#”号，那么“#”号后面的内容为注释.

10.1.2向量

R中最简单运算是向量的赋值，最常用的方法是使用函数c（），它把若干个数值或字符串组合为一个向量，比如：

>x1<-c（1,2）

>x1

[1]12

>x2<-c（3,4）

>x2

[1]34

>x<-c（x1,x2）

>x

[1]1234

可以对向量进行加（＋）、减（－）、乘（*）、除（/）和乘方（^）运算，表示对向量的每一个元素进行相应的运算.例如：

>x<-c（1,4,6.25）

>y<-x*2+1

>y

[1]3.09.013.5

也可以用向量作为函数自变量，sqrt、log、exp、sin、cos、tan等函数都可以用向量作自变量，结果是对向量的每一个元素取相应的函数值，如：

>sqrt（x）

[1]1.02.02.5

函数min和max分别取向量的最小值和最大值，函数sum计算自变量向量的元素和，函数range返回包含两个值的向量，第一个值是最小值，第二个值是最大值.例如：

>max（x）

[1]6.25

sort（x）返回x的元素从小到大排序的结果向量.order（x）返回使得x从小到大排列的元素下标向量.例如：

>a<-c（1,2.3,0.4）

>sort（a）

[1]0.41.02.3

>order（a）

[1]312

在R中很容易产生一个等差数列.例如，1:

n产生从1到n的整数列，-2:

3产生从-2到3的整数列，5:

2产生反向的数列：

>1:

5

[1]12345

>-2:

3

[1]-2-10123

>5:

2

[1]5432

要注意1:

n-1不是代表1到n-1而是向量1:

n减去1，这是一个常犯的错误：

>n<-5

>1:

n-1

[1]01234

>1:

（n-1）

[1]1234

seq函数是更一般的等差数列函数.如果只指定一个自变量n>0，则seq（n）相当于1:

n.指定两个自变量时，第一自变量是开始值，第二自变量是结束值，如seq（-2,3）是从-2到3.R函数调用的一个很好的特点是它可以使用不同个数的自变量，自变量可以用“自变量名＝自变量值”的形式指定.例如，seq（-2,3）可以写成seq（from=-2,to=3）.可以用一个by参数指定等差数列的增加值，例如：

>seq（0,2,0.7）

[1]0.00.71.4

也可以写成seq（from=0,to=2,by=0.7）.有参数名的参数的次序可以颠倒，如：

>seq（0,by=0.7,to=2）

[1]0.00.71.4

也可以用length参数指定数列长度，如seq（from=10,length=5）产生10到14的向量.

另一个类似的函数是rep，它可以重复第一个自变量若干次，例如：

>rep（x,3）

[1]1.004.006.251.004.006.251.004.006.25

第一个参数名为x，第二个参数名为times（重复次数）.

R提供了十分灵活的访问向量元素和向量子集的功能.某一个元素只要用x[i]的格式访问，其中x是一个向量名，或一个取向量值的表达式，如：

>x

[1]1.004.006.25

>x[2]

[1]4

可以单独改变一个元素的值，例如：

>x[2]<-125

>x

[1]1.00125.006.25

事实上，R提供了四种方法来访问向量的一部分，格式为x[v]，x为向量或向量值的表达式，v是一个向量，取值在1到length（x）之间，取值允许重复，例如：

>x[c（1,3）]

[1]1.006.25

>x[1:

2]

[1]1125

>x[c（1,3,2,1）]

[1]1.006.25125.001.00

向量v的取值在－length（x）到－1之间，表示扣除相应位置的元素.例如：

>x[-（1:

2）]

[1]6.25

10.1.3数组与矩阵

数组（array）可以看成是带多个下标的元素的集合，常用的是数值型的数组如矩阵，也可以有其它类型（如字符型、逻辑型、复型数组）.R可以很容易地生成和处理数组，特别是矩阵（二维数组）.

数组有一个特征属性叫做维数向量（dim属性），维数向量是一个元素取正整数值的向量，其长度是数组的维数，比如维数向量有两个元素时数组为二维数组（矩阵）.

一组值只有定义了维数向量（dim属性）后才能被看作是数组.比如：

>z<-1:

12

>dim（z）<-c（3,4）

>z

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

[1,]14710

[2,]25811

[3,]36912

这时z已经成为了一个维数向量为c（3,4）的二维数组.也可以把向量定义为一维数组，例如：

>dim（z）<-12

>z

[1]123456789101112

数组元素的排列次序缺省情况下按列排序，即第一下标变化最快，最后下标变化最慢，对于矩阵（二维数组）则是按列存放.例如，假设数组a的元素为1:

24，维数向量为c（2,3,4），则各元素次序为a[1,1,1],a[2,1,1],a[1,2,1],a[2,2,1],a[1,3,1],...,a[2,3,4].

用函数array（）或matrix（）可以更直观地定义数组.如

>array（1:

12,dim=c（3,4））

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

[1,]14710

[2,]25811

[3,]36912

函数matrix（）用来定义最常用的二维数组，即矩阵.例如，定义一个3行4列，由1:

12按行次序排列的矩阵，可以用：

>matrix（1:

12,ncol=4,byrow=T）

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

[1,]1234

[2,]5678

[3,]9101112

注意在有数据的情况下只需指定行数或列数之一.指定的数据个数允许少于所需的数据个数，这时循环使用提供的数据.例如：

>b<-matrix（0,nrow=3,ncol=4）

>b

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

[1,]0000

[2,]0000

[3,]0000

要访问数组的某个元素，只要象上面那样写出数组名和方括号内用逗号分开的下标即可，如a[2,1,2]，b[2,3].

数组可以进行四则运算（+，－,*,/，^），表示数组对应元素的四则运算，参加运算的数组一般应该是相同形状的（dim属性完全相同）.例如，假设A,B,C是三个形状相同的数组，则

>D<-C+2*A/B

计算得到的结果是A的每一个元素除以B的对应元素加上C的对应元素乘以2得到相同形状的数组.四则运算遵循通常的优先级规则.

矩阵是二维数组，但因为其应用广泛所以对它定义了一些特殊的运算和操作.函数t（A）返回矩阵A的转置.nrow（A）为矩阵A的行数，ncol（A）为矩阵A的列数.矩阵之间进行普通的加减乘除四则运算仍遵从一般的数组四则运算规则，即数组的对应元素之间进行运算，所以注意A*B不是矩阵乘法而是矩阵对应元素相乘.

要进行矩阵乘法，使用运算符%*%，A%*%B表示矩阵A乘以矩阵B（当然要求A的列数等于B的行数）.例如：

>A<-matrix（1:

12,nrow=4,ncol=3,byrow=T）

>B<-matrix（c（1,0）,nrow=3,ncol=2,byrow=T）

>A

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

[1,]123

[2,]456

[3,]789

[4,]101112

>B

[,1][,2]

[1,]10

[2,]10

[3,]10

>A%*%B

[,1][,2]

[1,]60

[2,]150

[3,]240

[4,]330

10.1.4数据框

数据框是R软件的一种数据结构.它通常是矩阵形式的数据，但矩阵各列可以是不同类型的.数据框每列是一个变量，每行是一个观测.例如

>air<-data.frame（data=c（1.13,1.16,1.18,1.19）,

+condition=c（"Fog","Clear","Cloud","Fog"）,

+Oxon=c（10.3,NA,36.4,13.6））

>air

dataconditionOxon

11.13Fog10.3

21.16ClearNA

31.18Cloud36.4

41.19Fog13.6

一个矩阵可以用data.frame（）转换为一个数据框，如果它原来有列名，则其列名被作为数据框的变量名，否则系统自动为矩阵的各列起一个变量名.如

>A<-matrix（1:

12,nrow=4,ncol=3,byrow=T）

>data.frame（A）

X1X2X3

1123

2456

3789

4101112

10.1.5程序控制结构

R是一个表达式语言，其任何一个语句都可以看成是一个表达式.表达式之间以分号分隔或用换行分隔.表达式可以续行，只要前一行不是完整表达式（比如末尾是加减乘除等运算符，或有未配对的括号）则下一行为上一行的继续.

若干个表达式可以放在一起组成一个复合表达式，作为一个表达式使用.组合用大括号表示，如：

>{

>x<-15

>x

>}

R语言也提供了其它高级程序语言共有的分支、循环等程序控制结构.分支结构包括if结构：

if（条件）表达式1

或

if（条件）表达式1else表达式2

其中的“条件”为一个标量的真或假值，表达式可以是用大括号包围的复合表达式.有else子句时一般写成：

if（条件）{

表达式组……

}else{

表达式组………

}

这样的写法可以使else不至于脱离前面的if.注意“&&”表示“与”，它是一个短路运算符，即第一个条件为假时就不计算第二个条件.在条件中也可以用“||”（两个连续的竖线符号）表示“或”，它也是短路运算符，当第一个条件为真时就不再计算第二个条件.

循环结构中常用的是for循环，是对一个向量或列表的逐次处理，格式为“for（nameinvalues）表达式”，如若要计算1到5这五个数的和：

>s<-0;x<-1:

5

>for（iin1:

5）{

+s<-s+x[i]

+}

>s

[1]15

当然，如果只是要求各元素的和，只要调用sum（x）即可.从这里我们也可以看出，显式的循环经常是可以避免的，利用函数对每个元素计算值、使用sum等统计函数往往可以代替循环.因为循环在R中是很慢的，所以应尽可能避免使用显式循环.

10.1.6函数

R软件允许用户自己创建模型的目标函数，有的R函数存储为特殊的内部形式，并可以被进一步的调用.这样在使用时可以使语言更有利、更方便，而且程序也更美观.学习写自己的程序是用户学习使用R语言的主要方法之一.

事实上，R软件提供的绝大多数函数，如mean（）、var（）等，是系统编写人员写在R语言中的函数，与用户自己写的函数本质上没有多大差别.

函数定义的格式为

name<-function（arg_1,arg_2,…）expression

expression是R软件中的表达式，arg_1,arg_2,…表示函数的参数.调用函数的格式为name（expr_1,expr_2）.如

>f<-function（x,y）cos（x）/（1+y^2）

>f（0,1）

[1]0.5

§10.2用R软件求解概率论与数理统计问题

10.2.1常见分布的概率计算

在R软件中，提供了计算典型分布的分布函数、分布律或概率密度函数、分位数等各种函数.下面以正态分布为例，介绍各种函数的具体命令.

pnorm（x,mu,sigma）表示均值为mu，标准差为sigma的正态分布在x处的分布函数值；在x处的密度函数为dnorm（x,mu,sigma）；下侧alpha分位点为qnorm（alpha,mu,sigma），注意，R软件计算是下侧分位点，而不是上侧分位点，其上侧alpha分位点为qnorm（1-alpha,mu,sigma）；产生n个服从均值为mu，标准差为sigma的正态分布的命令为rnorm（n,mu,sigma）.

例1设

.

（1）求

；

（2）已知

，求

；

（3）已知

，求

.

解

（1）

，

，

R命令如下所示：

>pnorm（1,4,1.5）

[1]0.02275013

>pnorm（3,4,1.5）-pnorm（0,4,1.5）

[1]0.2486622

>1-pnorm（6,4,1.5）+pnorm（2,4,1.5）

[1]0.1824224

所以

，

，

（2）

，则

，

就是

的下侧0.5分位数，R命令为

>qnorm（0.5,4,1.5）

[1]4

所以

.

（3）因为

，所以

所以

，R中命令为

>qnorm（0.975）*1.5#省略mu和sigma两个参数表示标准正态分布

[1]2.939946

故

.

对于其他分布一样，在分布名称的前面加上p、d、q、r分别表示分布函数、密度函数、下侧分位数和生成随机数.具体分布名称如表10-1所示.

表10.1

分布

R软件中的名称

括号中附加参数

（0-1）分布b（1,p）

binom

1,p

二项分布b（n,p）

binom

n,p

泊松分布π（λ）

pois

λ

超几何分布h（n,M,N）

hyper

M,N-M,n

均匀分布U（a,b）

unif

a,b

指数分布Exp（λ）

exp

λ

正态分布N（μ,σ）

norm

μ,σ

卡方分布χ2（n）

chisq

n

t分布t（n）

t

n

F分布F（m,n）

f

m,n

例2某工厂生产的产品中废品率为0.005，任意取出1000件，计算

（1）其中至少有2件废品的概率；

（2）能以90%以上的概率保证废品件数不超过多少？

解

（1）令

表示1000件产品中的废品数，那么

，所求概率为

，R中命令为

>1-pbinom（1,1000,0.005）

[1]0.959909

所以其中至少有2件废品的概率为0.959909.

（2）依题意，求

，使得

，我们求

满足

即可.

>qbinom（0.9,1000,0.005）

[1]8

所以能以90%以上的概率保证废品件数不超过8.

10.2.2描述性统计量

R中计算样本均值、样本方差、样本标准差的命令分别是mean（x）、var（x）和sd（x），其中x是样本构成的向量，如

>x<-c（1,2,3,4,5）

>mean（x）

[1]3

>var（x）

[1]2.5

>sd（x）

[1]1.581139

计算k阶原点矩和k阶中心矩可以通过对向量x做变换得到，如

>mean（x^2）#样本的二阶原点矩

[1]11

>mean（（x-mean（x））^2）#样本的二阶中心矩

[1]2

10.2.3参数估计

1.矩估计

前面我们讲过，对未知参数进行估计，按问题性质的不同可以分为两类：

点估计和区间估计.按照构造统计量方法的不同，点估计又分为矩估计和最大似然估计.

R中没有专门用于矩估计的函数.实际上，只要能得到样本矩，就可以得到矩估计.计算样本的k阶原点矩或中心矩的程序是非常简单的.

设

，从前面第七章的内容可知，

和

的矩估计分别为

例如，首先生成100个服从

的正态分布，可以计算两个参数的矩估计，在R软件中的命令如下：

>x<-rnorm（100,2,3）

>mean（x）

[1]1.676401

>mean（（x-mean（x））^2）

[1]9.458816

得到

和

的矩估计分别为1.676401和9.458816，与实际值2和9很接近，事实上，随着样本量（即生成的随机数的个数）的增大，得到的

和

的矩估计值会越来越接近实际值.

2.最大似然估计

求最大似然估计需要把对数似然函数最大化，在R中需要利用求极小值的函数optimize（）或nlm（）.optimize（）的一般用法是

optimize（f=,interval=,lower=min（interval）,

upper=max（interval）,maximum=FALSE,…）

其中f是求极小的目标函数，interval是包含有极小的初始区间，lower是初始区间的左端点，upper是初始区间的右端点，maximum=FALSE表示求极小值，否则（maximum=TRUE）表示求函数的极大值.

例3设

是从参数为

的指数分布中抽出的样本，求

的最大似然估计.

解对数似然函数为

求

的最大似然估计需要求上面的对数似然函数的最大值，为此在R中先生成服从参数为2的指数分布，在R中的具体命令为

>n<-20

>lambda<-2

>x<-rexp（n,lambda）

>loglike<-function（lambda）n*log（lambda）-lambda*sum（x）

>optimize（loglike,c（0,2）,maximum=TRUE）

$maximum

[1]1.723709

$objective

[1]-9.110644

得到的

的最大似然估计值为1.723709，与

的实际值非常接近.

nlm（）函数的一般用法是nlm（f,p,…），其中f是求极小的目标函数，p是参数的初始值.由于函数nlm（）求的是极小值，所以应该求负对数似然函数的最小值.

例4设

是来自

的一个样本，求

和

的最大似然估计.

解对数似然函数为

首先在R中生成服从30个服从

的随机数，在R中的具体命令为

>n<-30

>x<-rnorm（n）

>nloglike<-function（p）

n/2*log（2*pi）+n/2*log（p[2]）+1/2/p[2]*sum（（x-p[1]）^2）

>nlm（nloglike,c（0.5,2））

$minimum

[1]45.33265

$estimate

[1]-0.24864601.2023759

$gradient

[1]-1.372058e-058.509664e-06

$code

[1]1

$iterations

[1]10

得到的

和

的最大似然估计值分别为-0.2486460和1.2023759，与

和

的实际值非常接近.

3.正态总体均值与方差的区间估计

在R中直接利用t.test和var.test检验函数对正态总体的参数进行区间估计.t.test函数可以估计一个正态总体的均值或两个正态总体均值差的置信区间，其具体格式为：

t.test（x,y=NULL,alternative=c（"two.sided","less","greater"）,

…,var.equal=FALSE,conf.level=0.95,...）

其中x,y表示样本，alternative表示置信区间的类型，取值"two.sided"、"less"、"greater"分别表示双侧置信区间、单侧置信上限和单侧置信下限.var.equal=FALSE表示两个样本的方差不等，var.equal=TURE表示两个样本的方差相等.conf.level是置信水平，默认值是0.95.

例5为估计一件物体的重量

，将其称了10次，得到的重量（单位：

kg）为：

10.1，10，9.8，10.5，9.7，10.1，9.9，10.2，10.3，9.9

假设所称出的物体重量服从

，求该物体重量

的置信水平为0.95的置信区间.

解这里求的是一个正态总体均值的区间估计，R命令为

>x<-c（10.1,10,9.8,10.5,9.7,10.1,9.9,10.2,10.3,9.9）

>t.test（x）

OneSamplet-test

data:

x

t=131.5854,df=9,p-value=4.296e-16

alternativehypothesis:

truemeanisnotequal