统计软件实践R操作结果报告.docx

统计软件实践R操作结果报告.docx

《统计软件实践R操作结果报告.docx》由会员分享，可在线阅读，更多相关《统计软件实践R操作结果报告.docx（138页珍藏版）》请在冰豆网上搜索。

统计软件实践R操作结果报告

统计软件实践课程作业

R语言与统计分析-汤银才

姓名：

XXX

班级：

14经管实验班

第三章概率与分布

§3.1随机抽样

从52张扑克牌中抽取4张对应的R命令为:

>sample（1:

52,4）

[1]4820724

抛一枚均匀的硬币10次在R中可表示为:

>sample（c（"H","T"）,10,replace=T）

[1]"H""H""H""T""T""T""T""T""T""H"

掷一棵骰子10次可表示为:

>sample（1:

6,10,replace=T）

[1]6421354444

一名外科医生做手术成功的概率为0.90,那么他做10次手术在R中可以表示为:

>sample（c（"成功","失败"）,10,replace=T,prob=c（0.9,0.1））

[1]"失败""失败""成功""成功""成功""成功""成功""成功""成功""成功"

若以1表示成功,0表示失败,则上述命令可变为:

>sample（c（1,0）,10,replace=T,prob=c（0.9,0.1））

[1]1111111111

§3.2排列组合与概率的计算

抽取的4张是有次序的,因此使用排列来求解.所求的事件（记为A）概率为：

>1/prod（52:

49）

[1]1.539077e-07

抽取的4张是没有次序的,因此使用组合数来求解.所求的事件（记为B）概率为：

>1/choose（52,4）

[1]3.693785e-06

§3.4R中内嵌的分布

1）查找分布的分位数,用于计算假设检验中分布的临界值或置信区间的置信限.例如，显著性水平为5%的正态分布的双侧临界值是：

>qnorm（0.025）

[1]-1.959964

>qnorm（0.975）

[1]1.959964

2）计算假设检验的p值.比如自由度df1的χ23.84时的χ2检验的p值为：

>1-pchisq（3.84,1）

[1]0.05004352

而容量为14的双边t检验的p值为：

>2*pt（-2.43,df=13）

[1]0.0303309

§3.5应用:

中心极限定理

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（N>100）{

+rug（sample（x,100））

+}

+else{

+rug（x）

+}

+}

+}

二项分布:

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）

泊松分布:

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）

均匀分布:

unif（0,1）

>op<-par（mfrow=c（2,2））

>limite.central（）

>par（op）

指数分布:

exp

（1）

>op<-par（mfrow=c（2,2））

>limite.central（rexp,distpar=1,m=1,s=1）

>par（op）

正态混合分布：

>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.3

3.3从正态分布N（100;100）中随机产生1000个随机数,

1）作出这1000个正态随机数的直方图;

2）从这1000个随机数中随机有放回地抽取500个,作出其直方图;

3）比较它们的样本均值与样本方差.

（1）

>x=rnorm（1000,100,100）

>hist（x,prob=T,main="normalmu=0,sigma=1"）

（2）

>y=sample（x,500,replace=T）

>hist（y,prob=T,main="normalmu=0,sigma=1"）

（3）

>mean（x）

[1]99.39372

>mean（y）

[1]110.4469

>var（x）

[1]10362.39

>var（y）

[1]9457.672

其中x（1000个随机数）较y（500个随机数）样本均值小，而样本方差较大。

第四章探索性数据分析

§4.1常用分布的概率函数图

二项分布

>n<-20

>p<-0.2

>k<-seq（0,n）

>plot（k,dbinom（k,n,p）,type='h',

+main='Binomialdistribution,n=20,p=0.2',xlab='k'）

泊松分布：

>lambda<-4.0

>k<-seq（0,20）

>plot（k,dpois（k,lambda）,type='h',

+main='Poissondistribution,lambda=5.5',xlab='k'）

几何分布：

>p<-0.5

>k<-seq（0,10）

>plot（k,dgeom（k,p）,type='h',

+main='Geometricdistribution,p=0.5',xlab='k'）

超级和分布

>N<-30

>M<-10

>n<-10

>k<-seq（0,10）

>plot（k,dhyper（k,N,M,n）,type='h',

+main='Hypergeometricdistribution,

+N=30,M=10,n=10',xlab='k'）

负二项分布：

>n<-10

>p<-0.5

>k<-seq（0,40）

>plot（k,dnbinom（k,n,p）,type='h',

+main='NegativeBinomialdistribution,

+n=10,p=0.5',xlab='k'）

正态分布：

>curve（dnorm（x,0,1）,xlim=c（-5,5）,ylim=c（0,.8）,

+col='red',lwd=2,lty=3）

>curve（dnorm（x,0,2）,add=T,col='blue',lwd=2,lty=2）

>curve（dnorm（x,0,1/2）,add=T,lwd=2,lty=1）

>title（main="Gaussiandistributions"）

>legend（par（'usr'）[2],par（'usr'）[4],xjust=1,

+c（'sigma=1','sigma=2','sigma=1/2'）,

+lwd=c（2,2,2）,

+lty=c（3,2,1）,

+col=c（'red','blue',par（"fg"）））

###############t分布###########

>curve（dt（x,1）,xlim=c（-3,3）,ylim=c（0,.4）,

+col='red',lwd=2,lty=1）

>curve（dt（x,2）,add=T,col='green',lwd=2,lty=2）

>curve（dt（x,10）,add=T,col='orange',lwd=2,lty=3）

>curve（dnorm（x）,add=T,lwd=3,lty=4）

>title（main="StudentTdistributions"）

>legend（par（'usr'）[2],par（'usr'）[4],xjust=1,

+c（'df=1','df=2','df=10','Gaussiandistribution'）,

+lwd=c（2,2,2,2）,

+lty=c（1,2,3,4）,

+col=c（'red','blue','green',par（"fg"）））

##################χ2分布############

>curve（dchisq（x,1）,xlim=c（0,10）,ylim=c（0,.6）,col='red',lwd=2）

>curve（dchisq（x,2）,add=T,col='green',lwd=2）

>curve（dchisq（x,3）,add=T,col='blue',lwd=2）

>curve（dchisq（x,5）,add=T,col='orange',lwd=2）

>abline（h=0,lty=3）

>abline（v=0,lty=3）

>title（main='ChisquareDistributions'）

>legend（par（'usr'）[2],par（'usr'）[4],xjust=1,

+c（'df=1','df=2','df=3','df=5'）,

+lwd=3,lty=1,

+col=c（'red','green','blue','orange'）

+）

F分布

>curve（df（x,1,1）,xlim=c（0,2）,ylim=c（0,.8）,lty=1）

>curve（df（x,3,1）,add=T,lwd=2,lty=2）

>curve（df（x,6,1）,add=T,lwd=2,lty=3）

>curve（df（x,3,3）,add=T,col='red',lwd=3,lty=4）

>curve（df（x,3,6）,add=T,col='blue',lwd=3,lty=5）

>title（main="Fisher'sF"）

>legend（par（'usr'）[2],par（'usr'）[4],xjust=1,

+c（'df=（1,1）','df=（3,1）','df=（6,1）',

+'df=（3,3）','df=（3,6）'）,

+lwd=c（1,2,2,3,3）,

+lty=c（1,2,3,4,5）,

+col=c（par（"fg"）,par（"fg"）,par（"fg"）,'red','blue'））

#############对数正态分布##############

>curve（dlnorm（x）,xlim=c（-.2,5）,ylim=c（0,1.0）,lwd=2）

>curve（dlnorm（x,0,3/2）,add=T,col='blue',lwd=2,lty=2）

>curve（dlnorm（x,0,1/2）,add=T,col='orange',lwd=2,lty=3）

>title（main="Lognormaldistributions"）

>legend（par（'usr'）[2],par（'usr'）[4],xjust=1,

+c（'sigma=1','sigma=2','sigma=1/2'）,

+lwd=c（2,2,2）,

+lty=c（1,2,3）,

+col=c（par（"fg"）,'blue','orange'））

柯西分布：

>curve（dcauchy（x）,xlim=c（-5,5）,ylim=c（0,.5）,lwd=3）

>curve（dnorm（x）,add=T,col='red',lty=2）

>legend（par（'usr'）[2],par（'usr'）[4],xjust=1,

+c（'Cauchydistribution','Gaussiandistribution'）,

+lwd=c（3,1）,

+lty=c（1,2）,

+col=c（par（"fg"）,'red'））

威布尔分布

>curve（dexp（x）,xlim=c（0,3）,ylim=c（0,2））

>curve（dweibull（x,1）,lty=3,lwd=3,add=T）

>curve（dweibull（x,2）,col='red',add=T）

>curve（dweibull（x,.8）,col='blue',add=T）

>title（main="WeibullProbabilityDistributionFunction"）

>legend（par（'usr'）[2],par（'usr'）[4],xjust=1,

+c（'Exponential','Weibull,shape=1',

+'Weibull,shape=2','Weibull,shape=.8'）,

+lwd=c（1,3,1,1）,

+lty=c（1,3,1,1）,

+col=c（par（"fg"）,par（"fg"）,'red','blue'））

珈码分布

>curve（dgamma（x,1,1）,xlim=c（0,5）,lwd=2,lty=1）

>curve（dgamma（x,2,1）,add=T,col='red',lwd=2,lty=2）

>curve（dgamma（x,3,1）,add=T,col='green',lwd=2,lty=3）

>curve（dgamma（x,4,1）,add=T,col='blue',lwd=2,lty=4）

>curve（dgamma（x,5,1）,add=T,col='orange',lwd=2,lty=5）

>title（main="Gammadistributions"）

>legend（par（'usr'）[2],par（'usr'）[4],xjust=1,

+c（'k=1（Exponentialdistribution）',

+'k=2','k=3','k=4','k=5'）,

+lwd=c（2,2,2,2,2）,

+lty=c（1,2,3,4,5）,

+col=c（par（'fg'）,'red','green','blue','orange'））

贝塔分布

>curve（dbeta（x,1,1）,xlim=c（0,1）,ylim=c（0,4））

>curve（dbeta（x,3,1）,add=T,col='green'）

>curve（dbeta（x,3,2）,add=T,lty=2,lwd=2）

>curve（dbeta（x,4,2）,add=T,lty=2,lwd=2,col='blue'）

>curve（dbeta（x,2,3）,add=T,lty=3,lwd=3,col='red'）

>curve（dbeta（x,4,3）,add=T,lty=3,lwd=3,col='orange'）

>title（main="Betadistributions"）

>legend（par（'usr'）[1],par（'usr'）[4],xjust=0,

+c（'（1,1）','（3,1）','（3,2）',

+'（4,2）','（2,3）','（4,3）'）,

+lwd=c（1,1,2,2,3,3）,

+lty=c（1,1,2,2,3,3）,

+col=c（par（'fg'）,'green',par（'fg'）,

+'blue','red','orange'））

§4.2直方图与密度函数的估计

例4.2.1从二项分布binom（100,0.9）中抽取容量为N=100000的样本.试作出它的直方图及核密度估计曲线。

>N<-100000

>n<-100

>p<-.9

>x<-rbinom（N,n,p）

>hist（x,

+xlim=c（min（x）,max（x））,probability=T,

+nclass=max（x）-min（x）+1,col='lightblue',

+main='Binomialdistribution,n=100,p=.5'）

>>lines（density（x,bw=1）,col='red',lwd=3）

例4.2.2从负二项分布nbinom（10,0.25）中抽取容量为N=100000的样本.试作出它的直方图及核密度估计曲线。

>x<-rnbinom（N,10,.25）

>hist（x,

+xlim=c（min（x）,max（x））,probability=T,

+nclass=max（x）-min（x）+1,col='lightblue',

+main='Negativebinomialdistribution,n=10,p=.25'）

>lines（density（x,bw=1）,col='red',lwd=3）

§4.3单组数据的描述性统计分析

直方图

>library（DAAG）

>data（possum）

>fpossum<-possum[possum$sex=="f",]

>par（mfrow=c（1,2））

>library（DAAG）

>data（possum）

>fpossum<-possum[possum$sex=="f",]

>par（mfrow=c（1,2））

>attach（fpossum）

>hist（totlngth,breaks=72.5+（0:

5）*5,

ylim=c（0,22）,xlab="totallength",

main="A:

Breaksat72.5,77.5⋯"）

>hist（totlngth,breaks=75+（0:

5）*5,

ylim=c（0,22）,xlab="totallength",

main="B:

Breaksat75,80⋯"）

茎叶图

>stem（fpossum$totlngth）

Thedecimalpointisatthe|

74|0

76|

78|

80|05

82|0500

84|05005

86|05505

88|0005500005555

90|5550055

92|000

94|05

96|5

框须图

>library（DAAG）

>data（possum）

>fpossum<-possum[possum$sex=="f",]

>boxplot（fpossum$totlngth）

正态性检验

>qqnorm（fpossum$totlngth,

main="NormalityCheckviaQQPlot"）

>qqline（fpossum$totlngth,col='red'）

>dens<-density（totlngth）

>xlim<-range（dens$x）;ylim<-range（dens$y）

>par（mfrow=c（1,2））

>hist（totlngth,breaks=72.5+（0:

5）*5,

xlim=xlim,ylim=ylim,

probability=T,xlab="totallength",

main="A:

Breaksat72.5,77.5..."）

>lines（dens,col=par（'fg'）,lty=2）

>m<-mean（totlngth）

>s<-sd（totlngth）

>curve（dnorm（x,m,s）,col='red',add=T）

>hist（totlngth,breaks=75+（0:

5）*5,

xlim=xlim,ylim=ylim,

probability=T,xlab="totallength",

main="B:

Breaksat75,80..."）

>lines（dens,col=par（'fg'）,lty=2）

>m<-mean（totlngth）

>s<-sd（totlngth）

>curve（dnorm（x,m,s）,col='red',add=T）

>x<-sort（totlngth）

>n<-length（x）

>y<-（1:

n）/n

>m<-mean（totlngth）

>s<-sd（totlngth）

>plot（x,y,type='s',main="empiricalcdfof"）

>curve（pnorm（x,m,s）,col='red',lwd=2,add=T）

总体描述

>summary（fpossum$totlngth）

Min.1stQu.MedianMean3rdQu.Max.

75.0085.2588.5087.9190.5096.50

>mean（fpossum$totlngth）

[1]87.90698

五数及样本分位数概括

>fivenum（fpossum