统计计算-R软件-课程设计

课程设计报告

课程设计名称统计计算

课程设计题目基于R软件的MCMC算法的实现

专业统计学

班级_____________________

学生姓名_____________________

学生学号___________________

指导教师___________________

2016年7月8日

课程设计目的

统计计算课程设计是在学习了《统计计算》课程之后，进行此课程设计，

是对这门课程的全面复习，也是这门课程理论知识的实践，是整个教学工作的

重要环节。

通过统计计算课程设计教学所要到达的目的是：以统计计算课程和理论知

识为根底，通过课程设计的实践，加强学生对所学相关课程的理解、掌握，训

练并提高R软件的使用、统计蒙特卡洛方法、独立解决问题的能力。

二.设计内容

1.运用R软件完成《习题一》L5题。

2.编制梯形求积公式和抛物线求积公式计算丁十公的程序。

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

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

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

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

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

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

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

三.程序及结果

L题目：运用R软件完成《习题一》L5题。

程序代码

xl<-c(23,20,18,29,43,35,32,40,29,26,24,26);/x2<-c(1000,1000,500,500,500,100,100,100,10

0,100,100,100,100,100,100);/x3<-c(70,70,70,70,10);

mean(xl);#均值

var(xl):力方差

sd(xl);#标准差

median(xl);#中位数

which(table(xl)==max(table(xl)));#众数

cv=sd(x1)/mean(x1);#变异系数

print(cv);

n=length(xl);

gl=sqrt(l/6*n)*sum(((xl-mean(xl))/sd(xl))"3);#标准偏度系数

print(gl);

g2=sqrt(n/24)*((1/n)*sum(((x1-mean(x1))/sd(x1))-4)-3);#标准峰度系数，sqrt为平方根函数

print(g2);

结果

(1)xl<-c(23,2(),18,29,43,35,32,40,29,26,24,26)

>mean(xl);#均值

[1]28.75

>var(xl);#方差

[1]58.38636

>sd(xl);#标准差

[1]7.641097

>median(xl);#中位数

[1]27.5

>which(table(xl)==max(table(xl)));#众数

2629

56

>cv=sd(xl)/mean(xl);#变异系数

>print(cv);

[1]0.2657773

>n=length(xl);

>gl=sqrt(l/6*n)*sum(((xl-mean(xl))/sd(xl))"3);#标准偏度系数

>print(gl);

[1]7.433059

>g2=sqrt(n/24)*((1/n)*sum(((xl-mean(x1))/sd(xl))"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(xl);#均值

[1]300

>var(xl);#方差

[1]107142.9

>sd(xl);#标准差

[1]327.3268

>median(xl);#中位数

[1]100

>which(table(xl)==max(table(xl)));#众数

100

1

>cv=sd(x1)/mean(x1);#变异系数

>print(cv);

[1]1.091089

>n=length(xl);

>gl=sqrt(1/6*n)*sum(((xl-mean(x1))/sd(x1))3);#标准偏度系数

>print(gl);

[1]28.4031

>g2=sqrt(n/24)*((1/n)*sum(((xl-mean(xl))/sd(xl))4)-3);#标准峰度系数

>print(g2);

[1]-0.07153775

>

(3)x3<-c(70,70,70,70,10)

>mean(xl);#均值

[1]58

>var(xl);#方差

[1]720

>sd(xl);#标准差

[1]26.83282

>median(xl);#中位数

[1]70

>which(table(xl)==max(table(xl)));#众数

70

2

>cv=sd(xl)/mean(xl);#变异系数

>print(cv);

[1]0.4626348

>n=length(xl);

>gl=sqrt(l/6*n)*sum(((xl-mean(xl))/sd(xl))3);#标准偏度系数

>print(gl);

[1]-4.898979

>g2=sqrt(n/24)*((1/n)*sum(((xl-mean(xl))/sd(xl))"4)-3);#标准峰度系数

>print(g2);

[1]-0.4199206

>

2.题目：编制梯形求积公式计算工捻公的程序。

程序代码

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=O;n=l;

Wl=tixing(m,n)

C=integrate(f=integral,lower=m,upper=n)#积分

print(C[l]$value-Wl)

结果

>print(C[l]$value-Wl)

[1]0.01157178

3.题目：编制抛物线求积公式计算工六公的程序。

程序代码：

o=function(x){

y=x/(4+x*2)

return(y)

)

simpson=function(a,b){#simpson求积公式

xl=(a+b)/2;

y1=o(a);y2=o(xl);y3=o(b);

h=(b-a)/2

y=h/3*(y1+4*y2+y3)

return(y)

}

a=0;b=l;

W2=simpson(a,b)

C=integrate(f=o,lower=a,upper=b)#积分

print(C[l]$value-W2)

结果：

>print(C[l]$value-W2)

[1]-0.0001929302

4.题目：编制牛顿法求解Inx+x-1=0程序。

程序代码：

niudun:function(x){

y=log(x)+x-l;

return(y)

)

dniudun=function(x){

y=l/x+l;

return(y)

1

m=100：

n=rep(l,m)#重复

for(iin1:(m-2)){

n'i+l]=n[i]-niudun(n[i])/dniudun(n[i])

i:(abs(n[i+1]-n[i])<1*exp(-5))(break}Sabs返回整形数据的绝对值

)

print(n[l:i])

结果：

>print(a[1:i])

[1]1

5.题目：编制二分法的求解lnx+x-1=0程序。

程序代码：

e:=function(x){

y=log(x)+x-l

return(y)

)

a=rep(0,100);b=rep(0,100);c=rep(0,100);

a'l]=0.2;b[l]=2;

for(iin1: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+l]=a[i];b[i+l]=c[i]}

if(fb*fc<0){a[i+l]=c[i];b[i+l]=b[i]}

i:(abs(a[i]-b[i])<l*exp(-5)){break;}

print(c(i,a[i+l],b[i+l],c[i]))

)

c'i+l]=(a[i+l]+b[i+l])/2;

print(c[i+lj)

结果：

[1]1.00.21.11.1

[1]2.000.651.100.65

[1]3.0000.8751.1000.875

[1]4.00000.98751.10000.9875

[1]5.000000.987501.043751.04375

[1]6.0000000.9875001.0156251.015625

[1]7.0000000.9875001.0015621.001562

[1]8.00000000.99453121.00156250.9945312

[1]9.00000000.99804691.00156250.9980469

>c[i+l]=(a[i+l]+b[i+l])/2;

>print(c[i+l])

[1]1.000684

6.题目：用直接抽样法产生[1,10]上的均匀随机数的随机数。

程序代码：

junyun=function(n,a,b){

R=runif(n,0,1);

x=rep(0,n)

for(iin1:n){

x[i]=a+(b-a)*R[i]

)

return(x);

)

n=1000;a=l;b=10;

y=junyun(n,a,b)

mean(y)

var(y)

sy=sort(y)

i=((l:n)-0.5)/n

u=qunif(i,1,10)

plot(sy,u)

abline(0,1)

结果：

>mean(y)

[1]5.422167

>var(y)

[1]6.579233

7.题目：用直接抽样法产生[1,10]上指数分布的随机数。

程序代码：

zhishu=function(n,b)(

R=runif(n,0,1)

x=rep(0,n)

for(iin1:n){

x；i]=-(l/b)*log(R[i])

)

return(x)

)

n=1000;b=5

y=zhishu(n,b)

mean(y)

var(y)

sy=sort(y)

i=((l:n)-0.5)/n

u=qexp(i,b)

plot(sy,u)

abline(0,1)

结果：

>mean(y)

[1]0.2089918

>var(y)

[1]0.04107592

>

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

程序代码：

erwei=function(n){

rl=runif(n,0,1)#生成均匀分布的随机数

r2=runif(n,0,1)#生成均匀分布的随机数

u=rep(0,n)

v=rep(0,n)

for(iin1:n){

u'i]=sqrt(-2*log(rl[i]))*cos(2*pi*r2[i])#用变换抽样法产生正态分布随机数

v.i]=sqrt(-2*log(rl[i]))*sin(2*pi*r2[i])#|］变换抽样法产生正态分布随机数

}

return(cbind(u,v))

)

n=1000;

y二(erwei(n))

u=y[,1]

mean(u)

var(u)

su=sort(u)

i=((l:n)-0.5)/n

u=qnorm(i,0,1)

plot(su,u)

结果：

[1]-0.03751145

>var(u)

[1]1.003879

9.题目：用GFI法产生gamma分布的随机数。

程序代码；

GFI=function(n,a){

z=rep(0,n);

for(iin1:n){

y=l;A=0;

while(y>A){

r=runif(l,0,1);

x=-log(r);

y=runif(l,0,1)

A=(x/exp(x+l))"(a-l);

if(y<=A){z[i]=a*x}

}

}

return(z)

)

n=1000;a=2;

x=GFI(n,a);

mean(x);

var(x);

z=sort(x);

i=((l:n)-0.5)/n

y=qgamma(i,2,1);

plot(z,y);

abline(0,1)

结果：

>mean(x);

[1]2.044267

>var(x);

[1]2.182517

10题目：用舍选抽样法产生beta分布的随机数。

程序代码：

shexuan=function(n,a,b){

z=rcp(0,n);

fordin1:n){

r=2;A=l;

while(r>A){

x=runif(l,0,1);

r=runif(l,0,1);

px=(gamma(a+b)/(gamma(a)*gamna(b)))*x~(a-1)*(1-x)(b-l>;

pxO=(gamma(a+b)/(gamma(a)*g&mma(b)))*((a~l)/(a+b-2))~(a-l)*((b-l)/(a+b-2))”(bT);

A=px/pxO;

if(r<=A){z[i]=x}

}

}

return(z);

)

n=1000;a=2;b=l;

x=shexuan(n,a,b);

mean(x);

var(x);

sx=sort(x)

i=((l:n)-0.5)/n

u=qbeta(i,a,b)

plot(sx,u)

结果：

>mean(x);

[1]0.6802473

>var(x);

[1]

0.05645551

11题目：用复合抽样法产生密度函数为/。)=。+2(1-小"€[0,1],〃=0.5的随

机数。

程序代码：

fuhc=function(n,a)(

r=runif(n,0,1)

u=runif(n,0,1)

v=runif(n,0,1)

w=runif(n,0,1)

x=rep(0,n)

for(iin1:n)(

if(r[i]<=a){x[i]=u[i]}

if(r[i]>a){x[i]=max(v[i],w[i])}

}

return(x