第一题：

参考P82页中的实验2程序解决P88页问题2.5，并给出实验2中程序改进的建议

1.编写函数，解决P88问题2.5

  • 假设检验问题如下：
    • H0 :信号出现顺序随机
    • H1 :信号出现顺序不随机
  • 不妨令置信水平α=0.05，进行双边的随机游程检验(runs test)：
runs_test<-function(x,alternative=c("two.sided","less","greater")){#编写随机游程检验函数
  if(!is.vector(x))#检查x是否是向量
    stop('x is not a vector')
  if(any(is.na(x)))
    stop('NAs in x')
  x_factor<-as.factor(x)#将向量x转化为因子
  if(length(levels(x_factor))!=2)#数据中应该只有两个不同类别
    stop('x does not contain dichotomous data')
  alternative<-match.arg(alternative)
  t1=levels(x_factor)[1]
  t2=levels(x_factor)[2]
  n=length(x_factor)#计算游程长度
  n0=sum(x_factor==t1)
  n1=sum(x_factor==t2)
  #计算游程个数
  r=1
  for(i in 2:n){
    if(x_factor[i]!=x_factor[i-1]){r=r+1}
  }
  #算游程个数的概率
  p=NULL
  for(i in 1:n){
    k=(i-i%%2)/2
    if(i%%2==0){p[i]=(2*choose(n1-1,k-1)*choose(n0-1,k-1))/choose(n,n1)
    }else{
      p[i]=(choose(n1-1,k-1)*choose(n0-1,k)+choose(n0-1,k-1)*choose(n1-1,k))/choose(n,n1)
    }
  }
  p1<-sum(p[1:r])#开始算p-value
  if(alternative=='two.sided')
    if(p1>0.5){
      p_value<-(1-p1)*2
    }else{
      p_value<-p1*2
    }
  else if(alternative=='less')
    p_value<-p1
  else if(alternative=='greater')
    p_value<-1-p1
  else stop("irregulat alternative")
  DNAME <-deparse(substitute(x))
  METHOD <- "Runs Test"
  structure(list(alternative=alternative,p.value=p_value,method=METHOD,data.name=DNAME),class="htest")
}
signal<-c(0,1,0,1,1,1,0,0,1,1,0,0,0,0,1,1,1,1,1,1,1,1,1,0,1,0,0,1,1,1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,1,0,1,1,0,0,1,1,1,0,1,0,1,0,0,0,1,0,0,1,0,1,0,1,0,0,0,0,0,0,0,0)
runs_test(signal)
## 
##  Runs Test
## 
## data:  signal
## p-value = 0.8004
## alternative hypothesis: two.sided

从随机游程检验中可以看出，由于p-value远远大于0.05，所以不能拒绝零假设（信号出现顺序随机）。

2.给出实验2中程序改进的建议

  • 实验2的程序用来展示n0和n1给定的情况下随机游程个数R的精确分布是可以的，然而如果用它来进行随机游程检验的话仍然需要修改，例如检验时我们并不需要输出“概率分布表”（特别是n1和n2大的时候），只需要摘取P(R≤r)的概率即可。
  • 实验2的for循环可以用向量化运算替换掉以提高程序运行效率。
  • 实验2的程序输入参数只需要n0和n1两个即可，n可以由n0和n1计算而得。
  • 实验2的实验思路和实验过程好像没有什么联系……

第二题：

改进t检验，使其可以用于任意一组数据的平均数检验，画出程序流程图

  • 普通的t检验程序如下：
#[用置信区间来做假设检验]给出参数真值的1-α置信区间
t_test_interval<-function(data,alpha=0.05){
  x_bar<-mean(data)
  x_sd<-sd(data)
  n<-length(data)
  m<-qt(1-alpha/2,n-1)*x_sd/sqrt(n)
  up_t<-x_bar+m#置信上界
  down_t<-x_bar-m#置信下界
  return(c(down_t,up_t))
}

  • 基于bootstrap改进的t检验：

    • bootstrap无需分布假设，在传统t检验的统计假设不满足的时候，bootstrap可以派上用场。
    • bootstrap即从初始样本重复随机替换抽样，生成一个或一系列待检验统计量的经验分布，无需假设一个特定的理论分布，即可生成统计量的置信区间，并能检验统计假设。
    • 当原始样本性质足够好，自助样本的经验分布收敛到总体分布，就可以利用自助样本对总体进行统计推断
    • 具体步骤：
      • 1.原来已经从总体中随机抽出n个样本，这n个样本具有代表性。已经算出样本均值和样本标准差
      • 2.从样本中有放回随机抽取n个“自助”样本。有些观测可能会被选中多次，而有的观测可能一直都不会选中。可算得自助样本的均值和自主样本的标准差
      • 3.计算并记录对应的自助t统计量
      • 4.重复2和3一千次，得到自助t统计量的经验分布(一般来说1000次，实际上可以增加重复次数)
      • 5.将1000个自助t统计量从小到大排序
      • 6.找出自助t统计量的2.5%和97.5%分位点和
      • 7.此时总体均值的95%置信区间为,若原假设μ=μ0在置信区间内，则接受原假设，否则拒绝原假设
    • 基于bootstrap（自助法）改进的t检验程序如下：
t_test_bootstrap<-function(data,times=1000,alpha=0.05){
  x_bar<-mean(data)
  x_sd<-sd(data)
  n<-length(data)
  M=matrix(rep(NA,times*n),nrow = times,ncol=n)
  for(i in 1:times){
    M[i,]<-sample(data,n,replace = T)
  }
  t_boot<-function(x){(mean(x)-x_bar)/sd(x)}
  t_vector<-apply(M,1,t_boot)
  #求t_vector的分位点
  up_t_b<-sort(t_vector)[times*(1-alpha/2)]
  down_t_b<-sort(t_vector)[times*alpha/2]
  up_t<-x_bar+up_t_b*x_sd#置信上界
  down_t<-x_bar+down_t_b*x_sd#置信下界
  return(c(down_t,up_t))
}
  • 基于bootstrap改进的t检验的程序流程图(蓝色阴影为所修改的环节)

重要性：基于bootstrap改进的t检验是用bootstrap得出的自助t统计量的分布来代替原来的t分布，以构造参数真值的1-α置信区间，进而得出假设检验的结果。当分布未知或者受到极端值的影响的时候，用bootstrap得出的自助t统计量的分布比t分布更适用于均值（方差未知）的假设检验。

  • 演示（以书例2.1楼盘均价为例）
#楼盘均价数据
data_1<-c(36,32,31,25,28,36,40,32,41,26,35,35,32,87,33,35)
#普通t检验——参数真值的置信区间
t_test_interval(data_1)
## [1] 28.95415 44.04585
#基于bootstrap改进的t检验——参数真值的置信区间
t_test_bootstrap(data_1)
## [1] 19.05742 41.74768

可见，尽管基于bootstrap改进的t检验所给出的置信区间的长度大于普通t检验所给出的置信区间的长度（说明基于bootstrap改进的t检验比较粗糙），但是基于bootstrap改进的t检验所给出的置信上界明显小于普通t检验所给出的置信上界，这说明基于bootstrap改进的t检验能够减少异常值的影响，而且由于基于bootstrap改进的t检验没有正态分布假设，适用范围更宽。

第三题：

环境监测数据

  • 数据准备
#导入数据
library(readxl)
air<-read_excel("E:\\文档\\102-大学\\003-大三\\非参数统计\\2016.10.19-第五次非参作业\\某城市质量数据.xlsx",col_names=T)
#计算pm2.5的S+,S-
pm2.5<-air$`PM2.5`
## Warning: Unknown column 'PM2.5'
s_plus_pm2.5<-sum(pm2.5>40.4)
s_minus_pm2.5<-sum(pm2.5<40.4)
#计算aqi的S+,S-
aqi<-air$`AQI`
## Warning: Unknown column 'AQI'
s_plus_aqi<-sum(aqi>100)
s_minus_aqi<-sum(aqi<100)
  • 单侧3/4分位数检验
#pm2.5的单侧3/4分位数检验
pbinom(s_minus_pm2.5,34,3/4)
## [1] 3.388132e-21
#aqi的单侧3/4分位数检验
pbinom(s_minus_aqi,34,3/4)
## [1] 3.388132e-21

检验结果说明，在显著性水平α=0.1下，我们能拒绝“40.4是pm2.5的3/4分位数”的原假设，然而不能拒绝“100是aqi的3/4分位数”的原假设。

  • 单侧中位数检验
#pm2.5的单侧中位数检验
pbinom(s_minus_pm2.5,34,1/2)
## [1] 5.820766e-11
#aqi的单侧中位数检验
pbinom(s_plus_aqi,34,3/4)
## [1] 3.388132e-21

检验结果说明，在显著性水平α=0.1下，我们不能拒绝“40.4是pm2.5的中位数”的原假设，然而能拒绝“100是aqi的中位数”的原假设。

  • 对比与区别：
    • 根据以上单侧3/4分位数检验和单侧中位数检验的结果，我们可以得出如下结论：
      • 我们有信心拒绝“40.4是pm2.5的3/4分位数”的假设，而我们不能拒绝“40.4是pm2.5的中位数”的假设。
      • 我们有信心拒绝“100是aqi的中位数”的假设，而我们不能拒绝“100是aqi的3/4分位数”的假设。