基于R语言的统计模拟——点估计

月遇从云0408

于 2024-10-01 17:26:54 发布

阅读量411

点赞数 17

分类专栏：统计模拟在R语言中的实现文章标签： r语言开发语言学习笔记

本文链接：https://blog.csdn.net/weixin_63990108/article/details/142672147

版权

统计模拟在R语言中的实现专栏收录该内容

3 篇文章 1 订阅

订阅专栏

在前文中，我们已经探讨了如何利用R语言对假设检验进行模拟评价的方法，但在统计推断的广阔领域中，点估计同样占据着举足轻重的地位。本节将深入探索如何借助模拟手段来有效评估参数估计的性能与效果。

点估计，作为统计推断的核心组成部分，其目标在于通过样本数据精准地估计总体参数的真实值。然而，估计的准确性和可靠性并非自然而然的结果，而是需要经过科学严谨的方法和步骤来验证。因此，本节将聚焦于利用模拟技术，对点估计的准确性和效率进行全面评估。

通过模拟，我们可以构建出与真实世界相似的数据场景，进而在这些模拟数据中实施点估计方法。随后，通过对比估计结果与预设的总体参数真实值，我们可以直观地评估出点估计方法的优劣。此外，模拟技术还允许我们探究不同样本量、不同估计方法以及不同数据分布对点估计效果的影响，从而为实际应用中的参数估计提供有力的理论支撑和实践指导。

一、点估计

（一）定义与原理

点估计（point estimation）是根据样本数据估计总体分布所含未知参数的真值，所得到的值称为估计值。由于样本统计量为数轴上某一点值，估计的结果也以一个点的数值表示，所以称为点估计。点估计的目标是找到一个数字（即估计值），使得这个数字最好地表示总体参数（如均值、方差等）的值。

（二）构建方法

点估计的构造方法多种多样，常见的有点估计的常用方法有矩估计法、顺序统计量法、最大似然法、最小二乘法等。这些方法的核心思想都是通过构建一个近似值的函数，使得在给定一个数据集中的一个点位置时，这个函数的值尽可能地靠近给定点位置的实际值。

（三）优良性准则

为了评价点估计的优劣，需要制定一定的优良性准则。这些准则通常包括无偏性、一致性、有效性等。其中，无偏性是指估计量的数学期望等于被估计的总体参数；一致性是指当样本量趋于无穷大时，估计量依概率收敛于总体参数；有效性则是指估计量的方差尽可能小，即估计结果更加稳定可靠。

因此，在统计模拟的框架下，我们通常会采纳无偏性（Bias）、有效性（方差）以及MSE这三项关键指标，来全面评估参数估计方法的性能表现。

Bias估计：

$\widehat{Bias}\hat{(\theta)}=(1/N)\sum_{j=1}^{N}\hat{\theta}(X^j)-\theta$

方差估计：

$Var(\hat{\theta})=(1/(N-1))\sum_{j=1}^N(\hat{\theta}(X^j)-\hat{\theta}^.))$

MSE：

$MSE(\hat{\theta})=(1/N)\sum_{j=1}^N(\hat{\theta}(X^j)-\theta)^2$

（四）与区间估计的关系

点估计与区间估计是总体参数估计问题的两种不同方法。点估计关注的是找到一个最佳的数值来估计总体参数，而区间估计则是关注找到一个区间，使得这个区间最有可能包含总体参数的真实值。在实际应用中，点估计和区间估计可以相互补充，以提供更全面的信息。

（五）注意事项

点估计的结果是一个具体的数值，因此它无法给出估计的误差范围。如果需要了解估计的误差范围，通常需要结合区间估计来进行。
点估计的优劣取决于所选用的估计方法和样本数据的质量。因此，在进行点估计时，需要选择合适的估计方法和确保样本数据的可靠性和代表性。

二、R语言的实现

eg：已知参数ρ∈[-1,1]且X,η~N（0,1），我们定义

$Y=\rho*X+\sqrt{1-\rho^2}\eta$

通过Monte Carlo模拟研究，用样本相关系数估计总体ρ的bias及MSE，其中

ρ=-0.8,-0.5,-0.2,0,0.2,0.5,0.8
样本量取10，50，100，200，50
模拟次数：10000

#设置种子数
set.seed(1234)

#参数设置
p <- c(-0.8,-0.5,-0.2,0,0.2,0.5,0.8)
num <- c(10,50,100,200,500)
Bias.cor <- list() #用于存储Bias结果
MSE.cor <- list() #用于存储MSE结果

for (i in 1:length(p)) {

  bias.cor <- numeric()
  mse.cor <- numeric()
  
  for (k in 1:length(num)) {
    
    cor <- numeric()
    
    iter = 10000 #模拟次数
    for (j in 1:iter) {
      X <- rnorm(num[k],0,1) #创建随机数
      n <- rnorm(num[k],0,1)
      cor[j] = cor(X,p[i]*X+sqrt(1-(p[i])^2)*n) #利用公式计算相关系数
    }
    
    bias.cor[k] <- abs(mean(cor) - p[i])  #计算不同样本量下的不同p值的Bias
    mse.cor[k] <- mean((cor-p[i])^2) #计算不同样本量下的不同p值的MSE
  }
  
  Bias.cor[[i]] <- bias.cor
  MSE.cor[[i]] <- mse.cor
}

#绘图
#Bias
par(mar=c(5,5,2,2),xaxs="i",yaxs="i")
plot(num,Bias.cor[[1]],lwd=4,type="l",bty="l",lty=1,xaxt="n",yaxt="n",xlab="",
     ylab="",col=1,ylim=c(0,0.006))
axis(1,las=1,cex.axis=1.3,lwd=1.6) 
axis(2,las=1,cex.axis=1.3,lwd=1.6)
for(i in 2:length(p)){
  lines(num,Bias.cor[[i]],col=i,lwd=4,lty=i)
}
title(xlab="样本量",font.lab=7,cex.lab=1.5,line=2.6)
title(main="不同参数Bias随样本量变化")
legend(400,0.006,c("p=-0.8","p=-0.5","p=-0.2","p=0","p=0.2","p=0.5",
                  "p=0.8"),lwd=4,col=c(1:7),lty=c(1:7),bty="n",xpd=TRUE,cex=1.1)

#MSE
par(mar=c(5,5,2,2),xaxs="i",yaxs="i")
plot(num,MSE.cor[[1]],lwd=4,type="l",bty="l",lty=1,xaxt="n",yaxt="n",xlab="",
     ylab="",col=1,ylim=c(0,0.01))
axis(1,las=1,cex.axis=1.3,lwd=1.6) 
axis(2,las=1,cex.axis=1.3,lwd=1.6)
for(i in 2:length(p)){
  lines(num,MSE.cor[[i]],col=i,lwd=4,lty=i)
}
title(xlab="样本量",font.lab=7,cex.lab=1.5,line=2.6)
title(main="不同参数MSE随样本量变化")
legend(400,0.01,c("p=-0.8","p=-0.5","p=-0.2","p=0","p=0.2","p=0.5",
                  "p=0.8"),lwd=4,col=c(1:7),lty=c(1:7),bty="n",xpd=TRUE,cex=1.1)