🍉CSDN小墨&晓末:https://blog.csdn.net/jd1813346972
个人介绍: 研一|统计学|干货分享
擅长Python、Matlab、R等主流编程软件
累计十余项国家级比赛奖项,参与研究经费10w、40w级横向
该篇文章以实例的形式演示了利用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 模拟实验二:似然估计量的渐进正态过程
背景:总体为:
p
(
x
,
θ
)
=
{
θ
x
θ
−
1
,
0
<
x
<
1
0
,
e
l
s
e
,
g
(
θ
)
=
1
θ
p(x,\theta)=\begin{cases}\theta x^{\theta-1},0<x<1\\ 0,else\end{cases},g(\theta)=\frac{1}{\theta}
p(x,θ)={θxθ−1,0<x<10,else,g(θ)=θ1
用以上的方法给出
g
(
θ
)
=
1
θ
g(\theta)=\frac{1}{\theta}
g(θ)=θ1的似然估计量的渐进正态过程。
运行程序:
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 模拟实验三:方差估计量的优劣性比较
背景:正态总体 N ( 0 , σ 2 ) N(0,\sigma ^2) N(0,σ2),对于方差 σ 2 \sigma ^2 σ2的三个估计量:
σ 1 2 ^ = s 2 = 1 n − 1 ∑ i = 1 n ( x i − x ˉ ) 2 \hat{\sigma_1 ^2}=s^2=\frac{1}{n-1}\sum_{i=1}^n(x_i-\bar{x})^2 σ12^=s2=n−11∑i=1n(xi−xˉ)2, σ 2 2 ^ = 1 n ∑ i = 1 n x i 2 \hat{\sigma_2 ^2}=\frac{1}{n}\sum_{i=1}^nx_i^2 σ22^=n1∑i=1nxi2, σ 3 2 ^ = 1 n + 2 ∑ i = 1 n x i 2 \hat{\sigma_3 ^2}=\frac{1}{n+2}\sum_{i=1}^nx_i^2 σ32^=n+21∑i=1nxi2
利用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)
结果:
由表格及结果图可知, σ 1 2 \sigma_1^2 σ12最接近方差,即为最优估计量。