统计计算R软件课程设计

1、我春X 4大冬课程设计报告课程设计名称统计计算课程设计题目 基于R软件的MCMC：法的实现专 业 统计学班 级学生姓名 学生学号指导教师2016年 7月8 日课程设计目的统计计算课程设计是在学习了统计计算课程之后，进行此课 程设计，是对这门课程的全面复习，也是这门课程理论知识的实践， 是整个教学工作的重要环节。通过统计计算课程设计教学所要达到的目的是： 以统计计算课程 和理论知识为基础，通过课程设计的实践，加强学生对所学相关课程 的理解、掌握，训练并提高R软件的使用、统计蒙特卡洛方法、独立 解决问题的能力。二.设计内容.运用R软件完成习题一1.5题。 .一 .1.编制梯形求积公式和抛物线求积公

2、式计算 1：dx的程序。0 4 x2.编制二分法和牛顿法的求解lnx + x _1=0程序。.用直接抽样法产生1,10上的均匀随机数和指数分布的随机 数。.用二维变换抽样法产生正态分布的随机数。.用GFI法产生gamm吩布的随机数。.用舍选抽样法产生beta分布的随机数。.用复合抽样法产生密度函数为f (x) = a +2(1 -a)x,x w 0,1, a = 0.5的随机数。.用随机投针试验方法求兀的近似值。三.程序及结果.题目：运用R软件完成习题一1.5题。程序代码x1-c(23,20,18,29,43,35,32,40,29,26,24,26);/x2-c(1000,1000,500,

3、500,500,100 ,100,100,100,100,100,100,100,100,100);/x3-c(70,70,70,70,10);mean(x1);# 均值var(x1);# 方差sd(x1);#标准差median(x1);# 中位数which(table(x1)=max(table(x1); # 众数cv=sd(x1)/mean(x1);#变异系数print(cv);n=length(x1);g1=sqrt(1/6*n)*sum(x1-mean(x1)/sd(x1)A3);#标准偏度系数print(g1);sqrt 为g2=sqrt(n/24)*(1/n)*sum(x1-mea

4、n(x1)/sd(x1)A4)-3);#标准峰度系数，平方根函数print(g2) ;结果x1 mean(x1);# 均值1 28.75 var(x1);# 方差1 58.38636 sd(x1);# 标准差1 7.641097 median(x1);# 中位数27.5 which(table(x1)=max(table(x1); # 众数26 295 6cv=sd(x1)/mean(x1);#变异系数print(cv);1 0.2657773n=length(x1);g1=sqrt(1/6*n)*sum(x1-mean(x1)/sd(x1)A3);#标准偏度系数print(g1);1 7.4

5、33059g2=sqrt(n/24)*(1/n)*sum(x1-mean(x1)/sd(x1)A4)-3);#标准峰度系数print(g2);1 -0.7415435(2) x2 mean(x1);# 均值1 300 var(x1);# 方差1 107142.9 sd(x1);# 标准差1 327.3268 median(x1);# 中位数1 100which(table(x1)=max(table(x1); # 众数 1001cv=sd(x1)/mean(x1);# 变异系数print(cv);1 1.091089n=length(x1);g1=sqrt(1/6*n)*sum(x1-mean

6、(x1)/sd(x1)A3);#标准偏度系数print(g1);28.4031g2=sqrt(n/24)*(1/n)*sum(x1-mean(x1)/sd(x1)A4)-3);#标准峰度系数print(g2) ;1 -0.07153775(3) x3 mean(x1);# 均值1 58 var(x1);# 方差1 720 sd(x1);# 标准差1 26.83282 median(x1);# 中位数1 70which(table(x1)=max(table(x1); # 众数 702cv=sd(x1)/mean(x1);# 变异系数print(cv);1 0.4626348n=length(x

7、1);g1=sqrt(1/6*n)*sum(x1-mean(x1)/sd(x1)A3);#标准偏度系数print(g1);1 -4.898979g2=sqrt(n/24)*(1/n)*sum(x1-mean(x1)/sd(x1)A4)-3);#标准 峰度系数print(g2);1 -0.419920612.题目：编制梯形求积公式计算jdx的程序。0 4 x2程序代码integral=function(x)y=x/(4+xA2)return(y)tixing=function(m,n) #梯形求积公式ym= integral (m);yn= integral (n);y=0.5*(n-m)*(y

8、m+yn);return(y)m=0;n=1;W1=tixing(m,n)C=integrate(f= integral,lower=m,upper=n) # 积分print(C1 $value-W1)结果 print(C1 $value-W1)1 0.011571783.题目：编制抛物线求积公式计算f/Tdx的程序。 04 x程序代码：o=function(x)y=x/(4+xA2)return(y)simpson=function(a,b) #simpson求积公式x1=(a+b)/2;y1=o(a);y2=o(x1);y3=o(b);h=(b-a)/2y=h/3*(y1+4*y2+y3)

9、return(y)a=0;b=1;W2=simpson(a,b)C=integrate(f=o,lower=a,upper=b) # 积分print(C1 $value-W2)结果： print(C1 $value-W2)-0.00019293024.题目：编制牛顿法求解lnx+x_1=0程序。程序代码：niudun=function(x)y=log(x)+x-1;return(y)dniudun=function(x)y=1/x+1;return(y)m=100;n=rep(1,m) # 重复for(i in 1:(m-2)ni+1=ni-niudun(ni)/dniudun(ni)if(a

10、bs(ni+1-ni) print(a1:i)1 15.题目：编制二分法的求解lnx+x_1 = 0程序。程序代码：ef=function(x)y=log(x)+x-1return(y)a=rep(0,100);b=rep(0,100);c=rep(0,100);a1=0.2;b1=2;for(i in 1:100)ci=(ai+bi)/2;fa=ef(ai);fb=ef(bi);fc=ef(ci);if(fa*fc0) ai+1=ai;bi+1=ciif(fb*fc0) ai+1=ci;bi+1=biif(abs(ai-bi) mean(y)5.422167 var(y)1 6.579233

11、u7.题目：用直接抽样法产生1,10上指数分布的随机数。程序代码：zhishu=function(n,b) R=runif(n,0,1) x=rep(0,n) for(i in 1:n) xi=-(1/b)*10g(Ri) return(x) n=1000;b=5 y=zhishu(n,b) mean(y) var(y) sy=sort(y) i=(1:n)-0.5)/n u=qexp(i,b) plot(sy,u) abline(0,1)结果： mean(y) 1 0.2089918 var(y)1 0.04107592sy8.题目：用二维变换抽样法产生正态分布的随机数。程序代码:erwei

12、=function(n) r1=runif(n,0,1) r2=runif(n,0,1) u=rep(0,n) v=rep(0,n)for(i in 1:n)生成均匀分布的随机数生成均匀分布的随机数ui=sqrt(-2*log(r1i)*cos(2*pi*r2i) vi=sqrt(-2*log(r1i)*sin(2*pi*r2i) return(cbind(u,v)n=1000;y=(erwei(n)u=y,1mean(u)var(u)su=sort(u)i=(1:n)-0.5)/nu=qnorm(i,0,1) plot(su,u)用变换抽样法产生正态分布随机数#用变换抽样法产生正态分布随机数

13、结果:1 -0.03751145 var(u)1 1.0038799.题目：用GFI法产生gamm吩布的随机数。程序代码：GFI=function(n,a)z=rep(0,n);for(i in 1:n)y=1;A=0;while(yA)r=runif(1,0,1);x=-log(r);y=runif(1,0,1)A=(x/exp(x+1)A(a-1);if(yA)x=runif(1,0,1);r=runif(1,0,1);px=(gamma(a+b)/(gamma(a)*gamma(b)*xA(a-1)*(1-x)A(b-1);px0=(gamma(a+b)/(gamma(a)*gamma(

14、b)*(a-1)/(a+b-2)A(a-1)*(b-1)/(a+b-2)A(b-1)；A=px/px0;if(r=A) zi=x return(z); n=1000;a=2;b=1; x=shexuan(n,a,b); mean(x);var(x);sx=sort(x) i=(1:n)-0.5)/n u=qbeta(i,a,b) plot(sx,u)结果：mean(x);1 0.6802473var(x);10.0564555111 题目：用复合抽样法产生密度函数为f (x) =a +2(1 -a)x, x 0,1, a =0.5 的随机数。程序代码：fuhe=function(n,a)r=r

15、unif(n,0,1)u=runif(n,0,1)v=runif(n,0,1)w=runif(n,0,1)x=rep(0,n)for(i in 1:n)if(ria) xi=max(vi,wi)return(x)n=1000;a=0.5;x= fuhe (n,a);mean(x);var(x);f=function(y)z=(-a+sqrt(aA2+4*(1-a)*y)/(2*(1-a) return(z)i=(1:n)-0.5)/n;sx=sort(x);u=f(i);plot(sx,u);abline(0,1)结果： mean(x);1 0.579573 var(x);0.0781204412.题目：用随机投针试验方法求兀的近似值程序代码：touzhen=function(n,a,l)x=runif(n,0,0.5*a)y=runif(n,0,pi)S=0for(i in 1:n)if (xi=(l/2)*sin(yi) S=S+1 )pi2=2*l*n/(a*S)re