非参数Bootstrap

一些准备

knitr::opts_chunk$set(tidy = TRUE, 
                      warning = FALSE,
                      message = FALSE)

setwd("C:\\Users\\213yi\\Desktop\\非参数统计\\3-28")
library(showtext) #载入库

作业一

区间估计

data=read.table("nerve.txt")
fun<-function(data,alpha,fwd){
  data<-as.numeric(data)
  tn<-quantile(data,fwd)
  B=1000
  m=20
  assemble<-NULL
  resolution=as.data.frame(matrix(nrow=3,ncol = 3))
  colnames(resolution)<-c("下限","上限","区间长度")
  rownames(resolution)<-c("正太","枢轴量","分位数")
  for (i in 1:B){
    samp<-sample(data,m,replace = T)
    assemble<-c(assemble,quantile(samp,fwd))
  }
  #正太置信区间
  resolution[1,1]<-tn-qnorm(1-alpha/2)*sd(assemble)
  resolution[1,2]<-tn+qnorm(1-alpha/2)*sd(assemble)
  resolution[1,3]<-resolution[1,2]-resolution[1,1]
  #枢轴量置信区间
  resolution[2,1]<-2*tn-quantile(assemble,1-alpha/2)
  resolution[2,2]<-2*tn-quantile(assemble,alpha/2)
  resolution[2,3]<-resolution[2,2]-resolution[2,1]
  #分位数置信区间
  resolution[3,1]<-quantile(assemble,alpha/2)
  resolution[3,2]<-quantile(assemble,1-alpha/2)
  resolution[3,3]<-resolution[3,2]-resolution[3,1]
  
  return(resolution)
}

#0.75分位点 95%的置信区间
fun(data,0.05,0.75)

#0.75分位点 90%的置信区间
fun(data,0.10,0.75)

#0.25分位点 95%的置信区间
fun(data,0.05,0.25)

#0.25分位点 90%的置信区间
fun(data,0.10,0.25)

检验

#对正态性进行检验 
asse<-NULL
for (i in 1:1000){
    samp<-sample(data,20,replace = T)
    asse<-c(asse,as.numeric(quantile(samp,0.25)))
  }
shapiro.test(asse)#拒绝了正态性的假设
ks.test(asse,rnorm(100000,mean = mean(asse),sd=sd(asse)))#拒绝了正态性的假设
#H0:u=u0
qqnorm(asse)
qqline(asse)

• 因此，在假设条件满足且置信区间小的前提下，选择枢轴量置信区间

作业二

参数估计

方法：抽出100个点,估计出$\mu$=$\bar{X}$,从$e^{\mu }$中抽1000次,每次抽20个点,每次都计算$e^{\bar{X}}$，存储,枢轴量方法计算出θ

resolution=as.data.frame(matrix(nrow=2,ncol = 3))
colnames(resolution)<-c("下限","上限","区间长度")
rownames(resolution)<-c("参数bootstrap","非参数bootstrap")
sourced=rnorm(100,mean = 5,sd = 1)
xb=mean(sourced)

assum<-NULL
for (i in 1:1000){
  samp<-rnorm(20,mean=xb,sd=1)
  assum<-c(assum,exp(mean(samp)))
}
resolution[1,1]=2*exp(xb)-quantile(assum,0.975)
resolution[1,2]=2*exp(xb)-quantile(assum,0.025)
resolution[1,3]=resolution[1,2]-resolution[1,1]
df1<-assum

抽出100个点,作为数据源,抽1000次,每次抽20个点,每次都计算$e^{\bar{X}}$,存储,枢轴量方法计算出θ

assum<-NULL
for (i in 1:1000){
  samp<-sample(sourced,20)
  assum<-c(assum,exp(mean(samp)))
}
resolution[2,1]=2*exp(xb)-quantile(assum,0.975)
resolution[2,2]=2*exp(xb)-quantile(assum,0.025)
resolution[2,3]=resolution[2,2]-resolution[2,1]
df2<-assum
resolution

• 非参数Bootstrap方法比参数Bootstrap方法的置信区间长度更短，更为精确

可视化

#分布图
hist(df1,prob=TRUE,xlab="θ",ylab="密度",col="deepskyblue",main="")
hist(df2,prob=TRUE,xlab="",ylab="",col="red2",density=30,main="",add=TRUE)
legend("topright",legend=c("参数Bootstrap估计","非参数Bootstrap估计"),ncol=1,inset=0.04,col=c("deepskyblue","red2"),density=c(200,50),fill=c("deepskyblue","red2"),cex=0.8) 
title("Bootstrap估计直方图")

观察得到两者为偏态分布，猜想对数正态，并进行检验

检验

qqnorm(df1,main = "QQ-原始数据-参数")
qqline(df1)
qqnorm(log(df1),main = "QQ-log-参数")
qqline(log(df1))


qqnorm(df2,main = "QQ-原始数据-非参数")
qqline(df2)
qqnorm(log(df2),main = "QQ-log-非参数")
qqline(log(df2))

• 取对数之后，为正态分布
• 因此，数据为对数正态分布

作业三-正太记分检验

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-N5UNOeio-1650896609327)(C:/Users/213yi/Desktop/3.png)]

# 用正态记分检验方法评价人均年消费酒量的中位数
NScoreTest <- function(x, med.val, alternative = 'two-sided') {
  n <- length(x)
  x.diff <- x - med.val
  x.absdiff <- abs(x.diff)
  x.rank <- rank(x.absdiff)
  x.sign <- sign(x.diff)
  # 计算符号正态记分 x.nscore
  x.nscore <- qnorm((1 + (x.rank) / (n + 1)) / 2, 0, 1) * x.sign
  df <- data.frame(value = x, absvalue = x.absdiff, rank = x.rank,
  sign = x.sign, nscore = x.nscore)
  index <- order(x.absdiff)
  df <- df[index, ]
  w <- sum(x.nscore)
  Tval <- w / sqrt(sum(x.nscore ^ 2))
  if (alternative == 'two-sided') {
    pval = 2 * (1 - pnorm(abs(Tval), 0, 1))
  } else if (alternative == 'greater') {
    pval = 1 - pnorm(Tval, 0, 1) } else {
    pval = pnorm(Tval, 0, 1) }
    list(df = df, w = w, Tval = Tval, pval = pval)
}
alcohol <- c(4.12, 5.81, 7.63, 9.74, 10.39, 11.92, 12.32, 12.89, 13.54, 14.45)
NScoreTest(alcohol, med.val = 8, alternative = 'two-sided')

• p-value值为0.05567649,因此在显著性水平0.1的情况下，拒绝原假设，拒绝中位数为8的假设