该篇文章以实例的形式演示了利用R语言完成关于统计量近似分布的随机模拟:不同分布的正态性检验(正态分布、泊松分布);似然估计量的渐进正态过程;方差估计量的优劣性比较。
1 模拟实验一:不同分布的正态性检验
背景:总体分别为均匀分布,指数分布,两点分布(二阶矩存在)等情况下,利用随机模拟的方法画出样本均值的直方图,并进行正态性检验,给出p值
要求:(1)至少选两个总体,一个离散型,一个连续型。
(2)报告体现需出样本容量不同时,样本均值的趋近于正态分布的过程。
函数定义:
fenbuhanshu<- function (r,distpar,e,s,n,N)#r-随机数函数,distpar-参数取值范围,e-均值,s-标准差,n-样本数,N-抽取样本重复次数 #定义函数 { for (i in n) { 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*e )/(sqrt(i)*s) p=shapiro.test(x)$p.value #正态性检验 hist(x,col='green',probability=T,main=paste("n=",i,",p=",round(p,4)),ylim=c(0,max(.4, density(x)$y))) lines(density(x), col='red', lwd=3)}}
(1)连续型$N \sim U(0,1) $。
运行程序:
op=par(mfrow=c(2,2)) fenbuhanshu(runif,c(0,1),0.5,1/sqrt(12),c(3,30,300,500),1000) #分别做3,30,300,500次实验 par(op)
运行结果:
(2)离散型$N \sim P(1) $
运行程序:
op=par(mfrow=c(2,2)) limite.central(rpois,1,1,1,c(3,30,300,500),1000) par(op)
运行结果:
2 模拟实验二:似然估计量的渐进正态过程
背景:总体为:
运行程序:
t=par(mfrow=c(2,2)) #两行两列画布 f=function(theta){ logL=m*log(theta)+(theta-1)*sum(log(x)) return(logL)} #返回联合密度函数的ln值 n=c(5,50,500,5000) #分别模拟5,50,500,1000次 for (m in n){ a=c() for (i in 1:1000){ y=runif(m,min=0,max=1) x=sqrt(y) a[i]=1/(optimize(f,c(0,100),maximum = TRUE)$maximum)} #求最大似然函数值 p=shapiro.test(a)$p.value #Shapiro-Wilk方法进行正态检验 hist(a,col='green',probability=T,main=paste("n=",m,",p=",round(p,4)),ylim=c(0,max(.4, density(a)$y))) lines(density(a), col='red', lwd=3)} par(t)
运行结果:
由图可知,似然估计量随着样本越来越大有着渐进正态的过程。
3 模拟实验三:方差估计量的优劣性比较
利用R随机模拟方法,比较三个估计量分别在样本容量不同的情况下的优劣。给出表格和图形。
运行程序:
library(ggplot2) s1=c() s2=c() s3=c() op=par(mfrow=c(3,1)) #绘制3行1列画布 for (i in c(5,50,100)){ #样本容量分别设置为5,50,100 for (j in 1:100){ x=rnorm(i, mean=0,sd=1) #随机选取均值为0,方差为1的样本 y1=sum((x-mean(x))^2)/(i-1) s1[j]=y1 y2=sum(x^2)/i s2[j]=y2 y3=sum(x^2)/(i+2) s3[j]=y3} print(mean(s1)) print(mean(s1)-1) print(mean(s2)) print(mean(s2)-1) print(mean(s3)) print(mean(s3)-1) plot(density(s1),ylim=c(0,3),main=paste("n=",i)) lines(density(s2),col='red') lines(density(s3),col='green')} #分别绘制3根密度曲线 par(op)
结果:
由表格及结果图可知, 最接近方差,即为最优估计量。