统计计算 R软件 课程设计Word文档下载推荐.docx
《统计计算 R软件 课程设计Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《统计计算 R软件 课程设计Word文档下载推荐.docx（17页珍藏版）》请在冰豆网上搜索。
程序代码
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）））;
2629
56
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）
[1]300
[1]107142.9
[1]327.3268
[1]100
100
1
[1]1.091089
[1]28.4031
[1]-0.07153775
（3）x3<
-c（70,70,70,70,10）
[1]58
[1]720
[1]26.83282
[1]70
70
2
[1]0.4626348
[1]-4.898979
[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）;
}
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.01157178
3.题目:
编制抛物线求积公式计算
程序代码:
o=function（x）{
simpson=function（a,b）{#simpson求积公式
x1=（a+b）/2;
y1=o（a）;
y2=o（x1）;
y3=o（b）;
h=（b-a）/2
y=h/3*（y1+4*y2+y3）
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.0001929302
4.题目：
编制牛顿法求解
程序代码：
niudun=function（x）{
y=log（x）+x-1;
dniudun=function（x）{
y=1/x+1;
return（y）
m=100;
n=rep（1,m）#重复
for（iin1:
（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:
[1]1
5.题目：
编制二分法的求解
ef=function（x）{
y=log（x）+x-1
a=rep（0,100）;
b=rep（0,100）;
c=rep（0,100）;
a[1]=0.2;
b[1]=2;
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.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+1]=（a[i+1]+b[i+1]）/2;
print（c[i+1]）
[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=1;
b=10;
y=junyun（n,a,b）
mean（y）
var（y）
sy=sort（y）
i=（（1:
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）
x[i]=-（1/b）*log（R[i]）
return（x）
b=5
y=zhishu（n,b）
u=qexp（i,b）
[1]0.2089918
[1]0.04107592
8.题目：
用二维变换抽样法产生正态分布的随机数。
erwei=function（n）{
r1=runif（n,0,1）#生成均匀分布的随机数
r2=runif（n,0,1）#生成均匀分布的随机数
u=rep（0,n）
v=rep（0,n）
u[i]=sqrt（-2*log（r1[i]））*cos（2*pi*r2[i]）#用变换抽样法产生正态分布随机数
v[i]=sqrt（-2*log（r1[i]））*sin（2*pi*r2[i]）#用变换抽样法产生正态分布随机数
return（cbind（u,v））
y=（erwei（n））
u=y[,1]
mean（u）
var（u）
su=sort（u）
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）;
y=1;
A=0;
while（y>
A）{
r=runif（1,0,1）;
x=-log（r）;
y=runif（1,0,1）
A=（x/exp（x+1））^（a-1）;
if（y<
=A）{z[i]=a*x}
}
return（z）
a=2;
x=GFI（n,a）;
mean（x）;
var（x）;
z=sort（x）;
y=qgamma（i,2,1）;
plot（z,y）;
结果:
mean（x）;
[1]2.044267
var（x）;
[1]2.182517
10题目：
用舍选抽样法产生beta分布的随机数。
shexuan=function（n,a,b）{
z=rep（0,n）;
for（iin1:
r=2;
A=1;
while（r>
x=runif（1,0,1）;
r=runif（1,0,1）;
px=（gamma（a+b）/（gamma（a）*gamma（b）））*x^（a-1）*（1-x）^（b-1）;
px0=（gamma（a+b）/（gamma（a）*gamma（b）））*（（a-1）/（a+b-2））^（a-1）*（（b-1）/（a+b-2））^（b-1）;
A=px/px0;
if（r<
=A）{z[i]=x}
return（z）;
}
x=shexuan（n,a,b）;
sx=sort（x）
u=qbeta（i,a,b）
plot（sx,u）
[1]0.6802473
[1]0.05645551
11题目：
用复合抽样法产生密度函数为
fuhe=function（n,a）{
r=runif（n,0,1）
u=runif（n,0,1）
v=runif（n,0,1）
w=runif（n,0,1）
if（r[i]<
=a）{x[i]=u[i]}
if（r[i]>
a）{x[i]=max（v[i],w[i]）}
a=0.5;
x=fuhe（n,a）;
f=function（y）{
z=（-a+sqrt（a^2+4*（1-a）*y））/（2*（1-a））
n）-0.5）/n;
sx=sort（x）;
u=f（i）;
plot（sx,u）;
[1]0.579573
[1]0.07812044
12.题目：
用随机投针试验方法求π的近似值。
touzhen=function（n,a,l）{
x=runif（n,0,0.5*a）
y=runif（n,0,pi）
S=0
for（iin1:
if（x[i]<
=（l/2）*sin（y[i]））{S=S+1}
pi2=2*l*n/（a*S）
return（pi2）
n=1000000;
a=5;
l=2;
z=touzhen（n,a,l）
print（z）
[1]3.150152
四．参考资料
[1]《R统计建模与R软件》，薛毅，陈立萍.清华大学出版社，2007.
[2]《统计计算》，高惠璇编著，北京大学出版社，1995。
[3]《数值计算方法》，李乃成、邓建中编著，西安交通大学出版社，2002。
[4]《计算统计》，（美）吉文斯，（美）霍特伊著，王兆军，刘民千，邹长亮，杨建峰译，人民邮电出版社，2009。
[5]《统计模拟》，（美）罗斯著，王兆军，陈广雷，邹长亮译，人民邮电出版社，2007。
升级会员 