第五章 R数学计算与统计模拟

第五章R数学计算与统计模拟随着计算机技术的发展，统计模拟越来越得到重视。

R内置很多数学函数和统计分布函数，下面介绍这些函数及其应用以及如何编写简单的模拟实验。

一、随机抽样和数学计算（一）随机抽样在R中sample()函数可以用来产生随机排列或随机样本。

1.等可能不放回的随机抽样sample(x,n)其中，x为要抽取的向量，n为样本容量。

例如，从100个人中不放回地随机抽取10人的代码为：sample(1:100,10)如果要重复做抽样，并且抽样结果一样，可以设置随机种子：set.seed(10)设置种子的数字可以是任意的，只要固定、结果就一样。

2.等可能有放回的随机抽样sample(x,n,replace=TURE)其中选项replace=TURE表示抽样是有放回的，如果省略或为replace=FALSE表示抽样是不放回的。

例如，抛一枚硬币10次可以表示为：r<-sample(c("正","反"),10,replace = T)掷一颗骰子10次可以表示为：r<-sample(6,10,replace = T)注意，上面1:6简写为6。

3.不等可能的随机抽样sample(x,n,replace=TURE，prob=y)其中，prob=y用于指定x中元素出现的概率，向量y与x等长度。

例如，一名外科医生做手术成功的概率为0.9，那么他做10次手术成功的可能样本为：r<-sample(c("成功","失败"),10,replace = T,prob=c(0.9,0.1))R除了利用sample()函数产生进行随机抽样外，还可以利用随机分布函数产生多种不同分布下的随机数序列（见第二章）。

这些分布函数的形式为rfunc(n,p1,p2,…)，其中func指概率分布函数，n为生成数据的个数，p1, p2,…是分布的参数数值。

例如，模拟生成1000个自由度为2的卡方分布随机数，再求其平均数：mean(rchisq(1000,df=2))（二）数学计算1.一般数学计算一般数学计算函数参见第二章。

1)计算概率求向量所有元素的和：sum()，求向量所有元素的积：prod()。

例如，从一副完全打乱的52张扑克中抽取4张，若不考虑顺序，则4张全为A的概率为：p<-1/choose(52,4)若考虑顺序，则4张依次为红心A、方块A、黑桃A和梅花A的概率为：p<-1/prod(52:49)求向量所有元素的累计和：cumsum()，返回一个等长向量。

求向量所有元素的累计积：cumprod()，返回一个等长向量。

2)计算最（极）值求参数中所有元素的最小值：min()，求参数中所有元素的最大值：max()。

z<-matrix(1:6)min(z[,1],z[,2])求参数中所有向量对应元素的最小值：pmin()，返回一个向量；求参数中所有向量对应元素的最大值：pmax()，返回一个向量。

z<-matrix(1:6)pmin(z[,1],z[,2])以上函数是求向量的最值，如果要求函数的最值需要使用nlm()、optimize()和optim()等函数。

例如，求函数f(x)=x2-sin(x)的最小值，如下所示：fx=function(x){x^2-sin(x)}nlm(fx,3)函数nlm()求极小值，这里用到了数值分析的一种近似求根方法——Newton-type方法，第二个参数是用来设定初始值的。

也可以使用optimize()对一元函数最优化：optimize(fx,c(0,1))上面，c(0,1)是设定求最优化的区间。

默认求最小值，若要求最大值可以设置maximum=T。

多元函数的最优化问题比一元函数的复杂。

例如，求函数f(x,y)=100(y-x2)2+(1-x)2的极小值，如下所示。

绘制三维函数图：fxy=function(x,y){100*(y-x^2)^2+(1-x)^2}x=seq(-10,20,0.5)y=xz=outer(x,y,fxy)persp(x,y,z,theta=-30,phi=30,ticktype="detailed")其中，x,y都是长度为61的数值向量，outer(x,y,fxy)传回一个61×61的数值矩阵。

利用optim()函数求最优化，需要重新设定函数，原因是optim()函数要求输入函数的自变量是向量而非标量。

fx=function(x){100*(x[1]-x[2]^2)^2+(1-x[2])^2}optim(c(-10,-10),fx)函数预设求极小值，添加参数control=list(fnscale=-1)进行极大化。

线性规划问题可以用包lpSolve中的函数lp()求解，如求下面问题的最优化：#target: max C = 5×x1 + 8×x2#subject to:#x1 + x2 <= 2#x1 + 2×x2 = 3#x1,x2 >=0R语言代码如下：library(lpSolve)eg.lp <- lp(objective.in=c(5, 8),const.mat=matrix(c(1, 1, 1, 2), nrow=2),const.rhs=c(2, 3),const.dir=c("<=", "="), direction="max")eg.lp$solutioneg.lp$objval3)非线性方程求解R提供了一元函数求根的函数uniroot()：f=function(x){1-(1+x)^(-10)-6.71*x}plot(f,xlim=c(0.05,0.1))abline(h=0)x=uniroot(f,c(0.0001,1))x求解也可以化为最优化问题：g=function(x){f(x)^2}plot(g,xlim=c(0.05,0.1))optimize(g,c(0.0001,1))4)微积分R也可以进行微积分运算，包括符号微分和数值微分。

#求偏导数x <- 3; y <- 2.5; z <- 1exp1 <- expression(x / (y + exp(z)))exp1eval(exp1)D(exp1, "x")D(exp1, "y")D(exp1, "z")#求积分f<-function(x){x^2}integrate(f,0,1)另外，R的odesolve包也可以处理微分方程。

2.排序与配对1)排序对向量进行普通排序可以用sort()函数完成。

x<-c(1,5,3,2,8)sort(x)x注意，x本身没有改变。

要升序可加参数decreasing = FALSE。

如果想得到排序后的值在原向量中的索引，可以使用order()函数。

order(x)可以使用函数和索引对数据框进行排序。

v1<-c("f","m","f","f","m")#性别v2<-c(21,22,25,23,24)#年龄v<-cbind(v1,v2)r<-order(v2)rz<-v[r,]与排序相关的函数还有rank()，它返回向量中每一个元互的排位(rank)。

rank(x)注意rank()与order()的区别。

2)配对在d1和d2里有的是相同的记录，可以利用R的二元算子(binary operator)”%in%”找出。

”x%in%y”返回逻辑向量，表示在x中的元素是否也出现在y中。

1:7%in%5:105:10%in%1:7注意上面两条语句是不一样的。

如查找d1中和d2里名字相同的记录：d1[d1$姓名%in% d2$姓名,]在记录d中有重复记录，可以利用duplicated()函数找出。

v1<-c("f","m","f","f","m")#性别v2<-c(21,22,25,23,22)#年龄v<-cbind(v1,v2)duplicated(v)#后面与前面比较，返回逻辑向量v[duplicated(v),]3.向量和矩阵的计算1)一般运算向量乘以标量可以直接运算：x<-1:52*xR中有矩阵计算和处理的工具，函数rbind()和cbind()分别用上下或左右的方式合并向量或矩阵：m1 <- matrix(1, nr = 2, nc = 2)m2 <- matrix(2, nr = 2, nc = 2)rbind(m1, m2)数学意义上的两矩阵乘积的运算符是“%*%”，例如，对上面的矩阵m1和m2：m1%*%m2矩阵的转置由函数t()完成；这个函数也可以作用于数据框；函数det()可以用来求矩阵行列式的值；函数eigen()可以用来求矩阵的特征值和特征根；函数diag()可以用来提取或修正一个矩阵的对角元，或者创建一个对角矩阵。

函数sweep()可以做比较复杂的运算。

sweep(m1,1,c(1,2),"+")需要注意的是，diag()函数的再用性：如果它的参数是一个矩阵，它返回的是一个向量；反之亦然。

同样，如果它的参数是一个标量，那么这个函数返回指定大小的单位矩阵。

dm<-diag(m2)diag(dm)diag(3)x<-c(1,2,3,4,5)y<-c(2,4,6,8,10)2)内积和外积矩阵计算要区分内积和外积。

#内积也称点积x%*%y#或者crossprod(x,y)t(x)%*%y#外积、也称叉积x%o%y#或者outer(x,y)x%*%t(y)R中也有一些用于矩阵计算的专门的函数。

我们可以使用solve()求矩阵的逆，用qr()、chol()来分解矩阵，用svd()来做奇异值分解。

3)解线性方程组函数solve()还可以解线性方程组：#x1-x2=2#x1+x2=4a<-matrix(c(1,1,-1,1),2,2)b<-c(2,4)solve(a,b)4.集合运算R中包含一些方便的集合运算函数，主要有：集合x和y的并集：union(x,y)；集合x和y的交集：intersect(x,y)；集合x和y的差集：setdiff(x,y)，即所有属于集合x但不属于集合y的元素组成的集合；检验集合x和y是否相等：setequal(x,y)；检验元素c是否是集合x的元素：c%in%x；从含有n个元素的集合中选取含有k个元素的子集的数目：choose(n,k)。