R语言实验三

实验三数组的运算、求解方程（组）和函数极值、数值积分【实验类型】验证性【实验学时】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/yy^xexp(x)sqrt(y)x1-c(100,200);x2-1:6;x1+x22.x-1:5y-2*1:5x%*%ycrossprod(x,y)x%o%ytcrossprod(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)#提取对角线上的元素A[lower.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)x^3-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-1while(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(x[1]^2+x[2]^2-5,(x[1]+1)*x[2]-(3*x[1]+1))#定义函数组J-matrix(c(2*x[1],2*x[2],x[2]-3,x[1]+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)x^3-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*(x[2]-x[1]^2),1-x[1])#视为向量sum(F^2)}#向量对应分量平方后求和nlm(obj,c(-1.2,1))（2）fn-function(x){#定义目标函数F-c(10*(x[2]-x[1]^2),1-x[1])t(F)%*%F}#向量的内积gr-function(x){#定义梯度函数F-c(10*(x[2]-x[1]^2),1-x[1])J-matrix(c(-20*x[1],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-6N-1;h-b-a;T-h/2*(fun(a)+fun(b))#梯形面积repeat{h-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(-x^2)trape(f,-1,1)（3）常用求积分函数f-function(x)exp(-x^2)#定义函数integrate(f,0,1)integrate(f,0,10)integrate(f,0,100)integrate(f,0,10000)#当积分上限很大时，结果出现问题integrate(f,0,Inf)#积分上限为无穷大ft-function(t)exp(-(t/(1-t))^2)/(1-t)^2#对上述积分的被积函数e2作变量代换t=x/(1+x)后的函数integrate(ft,0,1)#与上述计算结果相同，且精度较高第二部分、教材例题：1.随机抽样（1）等可能的不放回的随机抽样: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(iinn){if(length(distpar)==2){x-matrix(r(i*N,distpar[1],distpar[2]),nc=i)}else{x-matrix(r(i*N,distpar),nc=i)}x-(apply(x,1,sum)-i*m)/(sqrt(i)*s)hist(x,col=lightblue,probability=T,main=paste(n=,i),ylim=c(0,max(.4,density(x)$y)))lines(density(x),col=red,lwd=3)curve(dnorm(x),col=blue,lwd=3,lty=3,add=T)if(N100){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-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)第三部分、课后习题：3.1a=sample(1:100,5)asum(a)3.2（1）抽到10、J、Q、K、A的事件记为A，概率为P(A)=1(5220)其中在R中计算得：1/choose(52,20)[1]7.936846e-15（2）抽到的是同花顺P(B)=(41)(91)(525)在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.33001]3.6op-par(mfrow=c(2,2))limite.central(rbeta,distpar=c(0.5,0.5),n=c(30,200,500,1000))par(op)