R语言实验三(共17页)

精选优质文档-倾情为你奉上实验三 数组的运算、求解方程（组）和函数极值、数值积分【实验类型】验证性【实验学时】2 学时【实验目的】1、掌握向量的四则运算和内积运算、矩阵的行列式和逆等相关运算；2、掌握线性和非线性方程（组）的求解方法，函数极值的求解方法；3、了解 R 中数值积分的求解方法。【实验内容】1、向量与矩阵的常见运算；2、求解线性和非线性方程（组）；3、求函数的极值，计算函数的积分。【实验方法或步骤】第一部分、课件例题：1. 向量的运算x-c(-1,0,2) y-c(3,8,2)v-2*x+y+1vx*yx/yyx exp(x)sqrt(y)x1-c(100,200); x2-1:6; x1+x22.x-1:5y-2*1:5x%*%y crossprod(x,y) x%o%y tcrossprod(x,y) outer(x,y)3.矩阵的运算A-matrix(1:9,nrow=3,byrow=T);AA+1 #A的每个元素都加上1B-matrix(1:9,nrow=3); BC-matrix(c(1,2,2,3,3,4,4,6,8),nrow=3); CD-2*C+A/B; D #对应元素进行四则运算x-1:9A+x #矩阵按列与向量相加E-A%*%B; E #矩阵的乘法y-1:3A%*%y #矩阵与向量相乘crossprod(A,B) #A的转置乘以Btcrossprod(A,B) #A乘以B的转置 4. 矩阵的运算A-matrix(c(1:8,0),nrow=3);At(A) #转置det(A) #求矩阵行列式的值diag(A) #提取对角线上的元素Alower.tri(A)=T-0;A #构造A对应的上三角矩阵qr.A-qr(A);qr.A #将矩阵A分解成正交阵Q与上三角阵R的乘积，该结果为一列表Q-qr.Q(qr.A);Q;R-qr.R(qr.A);R #显示分解后对应的正交阵Q与上三角阵Rdet(Q);det(R);Q%*%R #A=Q*Rqr.X(qr.A) #显示分解前的矩阵 5. 解线性方程组A-matrix(c(1:8,0),nrow=3,byrow=TRUE)b-c(1,1,1)x-solve(A,b); x #解线性方程组Ax=bB-solve(A); B #求矩阵A的逆矩阵BA%*%B #结果为单位阵6. 非线性方程求根f-function(x) x3-x-1 #建立函数uniroot(f,c(1,2) #输出列表中f.root为近似解处的函数值，iter为迭代次数，estim.prec为精度的估计值uniroot(f,lower=1,upper=2) #与上述结果相同polyroot(c(-1,-1,0,1) #专门用来求多项式的根，其中c(-1,-1,0,1)表示对应多项式从零次幂项到高次幂项的系数 7.求解非线性方程组（1）自编函数： (Newtons.R)Newtons-function (funs, x, ep=1e-5, it_max=100) index-0; k-1 while (k=it_max) #it_max 表示最大迭代次数 x1 - x; obj - funs(x); x - x - solve(obj$J, obj$f); #Newton 法的迭代公式 norm - sqrt(x-x1) %*% (x-x1) if (normep) index-1; break #index=1 表示求解成功 ; k-k+1 obj - funs(x); list(root=x, it=k, index=index, FunVal= obj$f) # 输出列表（2）调用求解非线性方程组的自编函数funs-function(x) f-c(x12+x22-5, (x1+1)*x2- (3*x1+1) # 定义函数组J-matrix(c(2*x1, 2*x2, x2-3, x1+1), nrow=2, byrow=T) # 函数组的 Jacobi 矩阵list(f=f, J=J) # 返回值为列表 : 函数值 f 和 Jacobi 矩阵 Jsource(F:/wenjian_daima/Newtons.R) # 调用求解非线性方程组的自编函数Newtons(funs, x=c(0,1)8.一元函数极值f-function(x) x3-2*x-5 # 定义函数optimize(f,lower=0,upper=2) # 返回值 : 极小值点和目标函数f-function(x,a) (x-a)2 # 定义含有参数的函数optimize(f,interval=c(0,1),a=1/3) # 在函数中输入附加参数9.多元函数极值（1）obj -function (x) # 定义函数 F-c(10*(x2-x12),1-x1) # 视为向量 sum (F2) # 向量对应分量平方后求和nlm(obj,c(-1.2,1)（2）fn-function(x) # 定义目标函数 F-c(10*(x2-x12), 1-x1) t(F)%*%F # 向量的内积gr - function(x) # 定义梯度函数 F-c(10*(x2-x12), 1-x1) J-matrix(c(-20*x1,10,-1,0),2,2,byrow=T) #Jacobi 矩阵 2*t(J)%*%F # 梯度optim(c(-1.2,1), fn, gr, method=BFGS) 最优点 (par) 、最优函数值 (value)10.梯形求积分公式（1）求积分程序： (trape.R)trape-function(fun, a, b, tol=1e-6) # 精度为 10 -6 N - 1; h - b-a ; T - h/2 * (fun(a) + fun(b) # 梯形面积 repeath - h/2; x-a+(2*1:N-1)*h; I -T/2 + h*sum(fun(x) if(abs(I-T) tol) break; N - 2 * N; T = I ; I（2）source(F:/wenjian_daima/trape.R) # 调用函数f-function(x) exp(-x2)trape(f,-1,1)（3）常用求积分函数f-function(x)exp(-x2) # 定义函数integrate(f,0,1)integrate(f,0,10)integrate(f,0,100)integrate(f,0,10000) # 当积分上限很大时，结果出现问题integrate(f,0,Inf) # 积分上限为无穷大ft sample(x, n) 其中x为要抽取的向量, n为样本容量（2）等可能的有放回的随机抽样: sample(x, n, replace=TRUE)其中选项replace=TRUE表示有放回的, 此选项省略或replace=FALSE表示抽样是不放回的sample(c(H, T), 10, replace=T)sample(1:6, 10, replace=T)（3）不等可能的随机抽样: sample(x, n, replace=TRUE, prob=y)其中选项prob=y用于指定x中元素出现的概率, 向量y与x等长度sample(c(成功, 失败), 10, replace=T, prob=c(0.9,0.1)sample(c(1,0), 10, replace=T, prob=c(0.9,0.1)2.排列组合与概率的计算1/prod(52:49)1/choose(52,4)3.概率分布qnorm(0.025) #显著性水平为5%的正态分布的双侧临界值qnorm(0.975)1 - pchisq(3.84, 1) #计算假设检验的p值2*pt(-2.43, df = 13) #容量为14的双边t检验的p值4. limite.central( )的定义limite.central - function (r=runif, distpar=c(0,1), m=.5, s=1/sqrt(12), n=c(1,3,10,30), N=1000) for (i in n) if (length(distpar)=2) x - matrix(r(i*N, distpar1,distpar2),nc=i) else x - matrix(r(i*N, distpar), nc=i) x 100 ) rug(sample(x,100) else rug(x) 5.直方图x=runif(100,min=0,max=1)hist(x)6.二项分布B（10，0.1）op - par(mfrow=c(2,2)limite.central(rbinom,distpar=c(10,0.1),m=1,s=0.9)par(op)7. 泊松分布: pios（1）op - par(mfrow=c(2,2)limite.central(rpois, distpar=1, m=1, s=1, n=c(3, 10, 30 ,50)par(op)8.均匀分布：unif(0,1)op - par(mfrow=c(2,2)limite.central( )par(op)9.指数分布：exp(1)op - par(mfrow=c(2,2)limite.central(rexp, distpar=1, m=1, s=1)par(op)10. 混合正态分布的渐近正态性mixn - function (n, a=-1, b=1)rnorm(n, sample(c(a,b),n,replace=T)limite.central(r=mixn, distpar=c(-3,3), m=0, s=sqrt(10), n=c(1,2,3,10)par(op)11. 混合正态分布的渐近正态性op - par(mfrow=c(2,2)mixn 1/choose(52,20)1 7.e-15（2）抽到的是同花顺P(B)= 在R中计算得： (choose(4,1)*choose(9,1)/choose(52,5)1 1.385e-053.3#(1)x-rnorm(1000,mean=100,sd=100)hist(x)#(2)y-sample(x,500)hist(y)#(3)mean(x)mean(y)var(x)var(y)3.4x-rnorm(1000,mean=0,sd=1)y=cumsum(x)plot(y,type = l)plot(y,type = p)3.5x-rnorm(100,mean=0,sd=1)qnorm(.025)qnorm(.975)t.test(x)由R结果知：理论值为-1.96,1.96，实际值为：-0.07929,0.330013.6 op - par(mfrow=c(2,2)limite.central(rbeta, distpar=c(0.5 ,0.5),n=c(30,200,500,1000)par(op)3.7N=seq(-4,4,length=1000)f-function(x) dnorm(x)/sum(dnorm(x)n=f(N)result=sample(n,replace=T,size = 1000)standdata=rnorm(1000)op-par(mfrow=c(1,2) #1行2列数组按列（mfcol）或行（mfrow）各自绘图hist(result,probability = T)lines(density(result),col=red,lwd=3)hist(standdata,probability = T)lines(density(standdata),col=red,lwd=3)par(op)专心-专注-专业