🍉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最接近方差，即为最优估计量。
