统计计算 R软件 课程设计

课程设计报告课程设计名称统计计算课程设计题目基于R软件的MCMC算法的实现专业统计学班级学生姓名学生学号指导教师2016年7 月8 日一． 课程设计目的统计计算课程设计是在学习了《统计计算》课程之后，进行此课程设计，是对这门课程的全面复习，也是这门课程理论知识的实践，是整个教学工作的重要环节。

通过统计计算课程设计教学所要达到的目的是：以统计计算课程和理论知识为基础，通过课程设计的实践，加强学生对所学相关课程的理解、掌握，训练并提高R 软件的使用、统计蒙特卡洛方法、独立解决问题的能力。

二． 设计内容1. 运用R 软件完成《习题一》1.5题。

2. 编制梯形求积公式和抛物线求积公式计算1204x dx x +⎰ 的程序。

3. 编制二分法和牛顿法的求解ln 10x x +-=程序。

4. 用直接抽样法产生[1,10]上的均匀随机数和指数分布的随机数。

5. 用二维变换抽样法产生正态分布的随机数。

6. 用GFI 法产生gamma 分布的随机数。

7. 用舍选抽样法产生beta 分布的随机数。

8. 用复合抽样法产生密度函数为()2(1),[0,1],0.5f x a a x x a =+-∈=的随机数。

9. 用随机投针试验方法求π的近似值。

三．程序及结果1.题目：运用R软件完成《习题一》1.5题。

程序代码x1<-c(23,20,18,29,43,35,32,40,29,26,24,26);/x2<-c(1000,1000,500,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))^3);#标准偏度系数print(g1);g2=sqrt(n/24)*((1/n)*sum(((x1-mean(x1))/sd(x1))^4)-3);#标准峰度系数，sqrt为平方根函数print(g2) ;结果（1）x1<-c(23,20,18,29,43,35,32,40,29,26,24,26)> mean(x1);#均值[1] 28.75> var(x1);#方差[1] 58.38636> sd(x1);#标准差[1] 7.641097> median(x1);#中位数[1] 27.5> which(table(x1)==max(table(x1))); #众数26 295 6> cv=sd(x1)/mean(x1);#变异系数> print(cv);[1] 0.2657773> n=length(x1);> g1=sqrt(1/6*n)*sum(((x1-mean(x1))/sd(x1))^3);#标准偏度系数> print(g1);[1] 7.433059> g2=sqrt(n/24)*((1/n)*sum(((x1-mean(x1))/sd(x1))^4)-3);#标准峰度系数> print(g2) ;[1] -0.7415435（2）x2<-c(1000,1000,500,500,500,100,100,100,100,100,100,100,100,100,100) > mean(x1);#均值[1] 300> var(x1);#方差[1] 107142.9> sd(x1);#标准差[1] 327.3268> median(x1);#中位数[1] 100> which(table(x1)==max(table(x1))); #众数1001> cv=sd(x1)/mean(x1);#变异系数> print(cv);[1] 1.091089> n=length(x1);> g1=sqrt(1/6*n)*sum(((x1-mean(x1))/sd(x1))^3);#标准偏度系数> print(g1);[1] 28.4031> g2=sqrt(n/24)*((1/n)*sum(((x1-mean(x1))/sd(x1))^4)-3);#标准峰度系数> print(g2) ;[1] -0.07153775>（3）x3<-c(70,70,70,70,10)> mean(x1);#均值[1] 58> var(x1);#方差[1] 720> sd(x1);#标准差[1] 26.83282> median(x1);#中位数[1] 70> which(table(x1)==max(table(x1))); #众数702> cv=sd(x1)/mean(x1);#变异系数> print(cv);[1] 0.4626348> n=length(x1);> g1=sqrt(1/6*n)*sum(((x1-mean(x1))/sd(x1))^3);#标准偏度系数> print(g1);[1] -4.898979> g2=sqrt(n/24)*((1/n)*sum(((x1-mean(x1))/sd(x1))^4)-3);#标准峰度系数> print(g2) ;[1] -0.4199206>2.题目：编制梯形求积公式计算1204x dx x +⎰ 的程序。

程序代码integral=function(x){y=x/(4+x^2)return(y)}tixing=function(m,n){ #梯形求积公式ym= integral (m);yn= integral (n);y=0.5*(n-m)*(ym+yn);return(y)}m=0;n=1;W1=tixing(m,n)C=integrate(f= integral,lower=m,upper=n) #积分print(C[1] $value-W1)结果> print(C[1] $value-W1)[1] 0.011571783.题目: 编制抛物线求积公式计算124x dx x +⎰ 的程序。

程序代码: o=function(x){y=x/(4+x^2)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)return(y)}a=0;b=1;W2=simpson(a,b)C=integrate(f=o,lower=a,upper=b) #积分print(C[1] $value-W2)结果：> print(C[1] $value-W2)[1] -0.00019293024.题目：编制牛顿法求解ln10x x+-=程序。

程序代码：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)){n[i+1]=n[i]-niudun(n[i])/dniudun(n[i])if(abs(n[i+1]-n[i])<1*exp(-5)) {break} #abs 返回整形数据的绝对值}print(n[1:i])结果：> print(a[1:i])[1] 15. 题目：编制二分法的求解ln10+-=程序。

x x程序代码：ef=function(x){y=log(x)+x-1return(y)}a=rep(0,100);b=rep(0,100);c=rep(0,100);a[1]=0.2;b[1]=2;for(i in 1:100){c[i]=(a[i]+b[i])/2;fa=ef(a[i]);fb=ef(b[i]);fc=ef(c[i]);if(fa*fc<0) {a[i+1]=a[i];b[i+1]=c[i]}if(fb*fc<0) {a[i+1]=c[i];b[i+1]=b[i]}if(abs(a[i]-b[i])<1*exp(-5)){break;}print(c(i,a[i+1],b[i+1],c[i]))}c[i+1]=(a[i+1]+b[i+1])/2;print(c[i+1])结果：[1] 1.0 0.2 1.1 1.1[1] 2.00 0.65 1.10 0.65[1] 3.000 0.875 1.100 0.875[1] 4.0000 0.9875 1.1000 0.9875[1] 5.00000 0.98750 1.04375 1.04375[1] 6.000000 0.987500 1.015625 1.015625[1] 7.000000 0.987500 1.001562 1.001562[1] 8.0000000 0.9945312 1.0015625 0.9945312[1] 9.0000000 0.9980469 1.0015625 0.9980469> c[i+1]=(a[i+1]+b[i+1])/2;> print(c[i+1])[1] 1.0006846. 题目：用直接抽样法产生[1,10]上的均匀随机数的随机数。

程序代码：junyun=function(n,a,b){R=runif(n,0,1);x=rep(0,n)for(i in 1:n){x[i]=a+(b-a)*R[i]}return(x);}n=1000;a=1;b=10;y=junyun(n,a,b)mean(y)var(y)sy=sort(y)i=((1:n)-0.5)/nu=qunif(i,1,10)plot(sy,u)abline(0,1)结果：> mean(y)[1] 5.422167> var(y)[1] 6.579233246810246810syu7. 题目：用直接抽样法产生[1,10]上指数分布的随机数。

程序代码：zhishu=function(n,b){R=runif(n,0,1)x=rep(0,n)for(i in 1:n){x[i]=-(1/b)*log(R[i])}return(x)}n=1000;b=5y=zhishu(n,b)mean(y)var(y)sy=sort(y)i=((1:n)-0.5)/nu=qexp(i,b)plot(sy,u)abline(0,1)结果：> mean(y)[1] 0.2089918> var(y)[1] 0.04107592>0.00.10.20.30.40.50.00.10.20.30.40.5syu8. 题目：用二维变换抽样法产生正态分布的随机数。