第一题：

画势函数图像,探究不同因素对势函数的影响

1.参数真值

  • 假设总体X来自Poisson分布P(λ),简单随机抽样X1,X2……Xn,假设检验问题H0:λ≥1,H1:λ<1.选取充分统计量ΣXi为检验统计量。我们在样本量为10时，对C=5和C=7考虑了检验势函数随分布 λ0从0变化到2的情况.在这里,拒绝域为{ΣXi<C}.
pois_power <- function(n,C,lty){ #n是样本量,C是拒绝域中的C,lty是画图线条类型
  lambda <- seq(from=0.1,to=2.0,by=0.01)
  pois_d <- ppois(C,n*lambda)
  plot(lambda,pois_d,type="l",xlab="λ",ylab="power",lty=lty,sub="检验势函数随分布参数变化曲线") #画势函数
  text(lambda[50],pois_d[50],paste("C=",as.character(C),sep = "")) #添加必要注释
}
pois_power(10,5,1) #画n=10,C=5时的势函数,线条类型是实线
par(new=T)
pois_power(10,7,2) #画n=10,C=7时的势函数,线条类型是虚线

2.样本大小[n=10和n=20]

pois_power_2 <- function(n,C,lty){ #n是样本量,C是拒绝域中的C,lty是画图线条类型
  lambda <- seq(from=0.1,to=2.0,by=0.01)
  pois_d <- ppois(C,n*lambda)
  plot(lambda,pois_d,type="l",xlab="λ",ylab="power",lty=lty,axes=F) #画势函数
  text(lambda[50],pois_d[50],paste("n=",as.character(n),sep = "")) #添加必要注释
}
pois_power_2(10,5,1) #画n=10,C=5时的势函数,线条类型是实线
par(new=T)
pois_power_2(20,5,2) #画n=20,C=5时的势函数,线条类型是虚线
axis(1)
axis(2,at=seq(0.2,1.0,0.2))
box()

3.显著性水平[α=0.1,0.05,0.01,0.005]

  • 假设总体X来自正态分布N(μ,σ2),简单随机抽样X1,X2……Xn,假设检验问题H0:μ=μ0,H1:μ>μ0.建立水平为α的检验.在原假设H0和标准差已知为σ0下,检验统计量.在这里,拒绝域为{u≥u1-α}.
normal_power<- function(miu0,sigma,alpha,n){ #miu0为μ0,sigma为σ0,alpha为显著性水平α,n为样本量
  miu<-seq(from=miu0-2*sigma,to=miu0+2*sigma,by=0.01)
  nor_d<- 1-pnorm(qnorm(1-alpha)-(miu-miu0)/(sigma/sqrt(n)))
  plot(miu,nor_d,type="l",xlab="μ",ylab="power")
  lines(c(miu0,miu0),c(0,alpha),lty=2)
  lines(c(0,miu0),c(alpha,alpha),lty=2)
  text(miu0-1.7*sigma,alpha*2,paste("α=",as.character(alpha),sep=""))
}

par(mfrow=c(2,2)) 
normal_power(3,1,0.1,25) #画μ0=3,σ0=1,显著性水平α=0.1,样本量为25时的势函数
normal_power(3,1,0.05,25) #画μ0=3,σ0=1,显著性水平α=0.05,样本量为25时的势函数
normal_power(3,1,0.01,25) #画μ0=3,σ0=1,显著性水平α=0.01,样本量为25时的势函数
normal_power(3,1,0.005,25) #画μ0=3,σ0=1,显著性水平α=0.005,样本量为25时的势函数

4.检验统计量的选择

  • 假设总体X来自指数分布Exp(λ),简单随机抽样X1,X2……Xn,假设检验问题H0:λ=λ0,H1:λ≠λ0.建立水平为α的检验.
  • 检验统计量的选择不止一种,在这里,我列出常见的两种检验统计量:
    • 2λ0Tn[其中Tn=ΣXi]
    • 2λ0nx(1)
  • 原假设H0成立时,上面两个检验统计量分别符合chisq(2n)和chisq(2),而原假设H0不成立时，上面两个检验统计量分别符合Gamma(n,λ/2λ0)和Gamma(1,λ/2λ0).
  • 两种检验统计量的势函数图像如下:
exp_pow_1 <- function(lambda0,n,alpha,lty){ #检验统计量为2λ0*Tn
  chisq1 <- qchisq(alpha/2,2*n)  #拒绝域左边界
  chisq2 <- qchisq(1-alpha/2,2*n)#拒绝域右边界
  lambda <- seq(lambda0*0.5,lambda0*1.5,by=0.01)
  power_exp <- pgamma(chisq1,shape=n,rate=lambda/(2*lambda0))+pgamma(chisq2,shape=n,rate=lambda/(2*lambda0),lower.tail=F)                     #λ不等于λ0时检验统计量落入拒绝域的概率
  plot(lambda,power_exp,type="l",lty=lty,xlim=c(1,3),ylim=c(0,0.6),xlab="λ")
  text(lambda[20],power_exp[20],"2λ0*Tn")
}
exp_pow_2 <- function(lambda0,n,alpha,lty){ #检验统计量为2λ0*n*x(1)
  chisq1 <- qchisq(alpha/2,2)    #拒绝域左边界
  chisq2 <- qchisq(1-alpha/2,2)  #拒绝域右边界
  lambda <- seq(lambda0*0.5,lambda0*1.5,by=0.01)
  power_exp <- pgamma(chisq1,shape=1,rate=lambda/(2*lambda0))+pgamma(chisq2,shape=1,rate=lambda/(2*lambda0),lower.tail=F)                     #λ不等于λ0时检验统计量落入拒绝域的概率
  plot(lambda,power_exp,type="l",lty=lty,xlim=c(1,3),ylim=c(0,0.6),xlab="λ")
  text(lambda[30],power_exp[30],"2λ0*n*x(1)")
}
exp_pow_1(2,10,0.05,1)          #画λ=2,n=10,α=0.05,检验统计量为2λ0*Tn时的势函数,线条类型是实线
par(new=T)  
exp_pow_2(2,10,0.05,2)          #画λ=2,n=10,α=0.05,检验统计量为2λ0*n*x(1)时的势函数,线条类型是虚线

  • 显然,用2λ0Tn的检验效率更高,这两种检验统计量之间应该选择2λ0Tn.

总结:我在第一题中分别画出泊松分布,正态分布和指数分布假设检验的势函数.影响势函数的四个因素:显著性水平,参数真值,样本大小和检验统计量的选择都已经通过画图比较的途径呈现在上方.

第二题：

(1)通读附件中的文章,叙述文章主要研究的问题,并归纳梳理文中传统的研究机理.

  • 文章主要研究的问题是:利用统计学的方法对某城市的土壤重金属污染状况进行分析研究,根据统计数据确定样本点的八种重金属元素(砷,镉,铜,镍,汞,铅,锌,铬)的浓度和样本点的地理位置来分析城市表层土壤重金属污染的主要原因.

  • 文中传统的研究机理:运用复杂的外部因素的机理分析方法来研究重金属污染,包括污染物的排放途径(采矿,工业“三废”排放,含重金属污水灌溉,城市生活垃圾,汽车尾气,污泥农用及含重金属的农药有机肥,使用重金属制品等人为因素)、污染物的传播途径(例如大气扩散,大气沉降,地表水扩散)、土壤或生物对污染物中重金属的吸附和释放机理等等.

(2)文章中的经验分布函数与教材中的有何不同?文章中这样构造的表2有何不足?表现在哪里?

  • 文章中的经验分布函数与教材中的有何不同?

    • 教材中的经验分布函数是,将容量为n的简单随机样本的值看作是某个离散随机变量等可能取的值而得到的离散分布的分布函数.(在这里这n个值的概率相等都为1/n)

    • 文章中的经验分布函数相当于将样本等距分割成若干组之后的累积分布函数,而每一组的概率由该组的频率决定.

    • 文章中构造经验分布函数的做法实际上损失了信息,将数值变量活生生地压缩成分类变量:原始数据能构造出文中的经验分布函数,文中的经验分布函数却推不出原始数据,而教材中的经验分布函数与原始数据是一一对应的.这样的恶果可由表2中体现出来.

  • 文章中这样构造的表2有何不足?表现在哪里?

    • 文章中说到,“虽然样本点的概率发生很小,然而他们会出现影响程度较大的浓度污染,这些浓度大的点对As污染土壤将回起到决定性因素”,然后说要找98%分位数及以上的数据来重点分析.文章说了这么多，实质上,就是挑出每一种重金属污染浓度最高的2%来进行重点分析而已.
    • 首先回答一下题目问题：实际上文中选取As浓度大于13.018以上的数据进行分析,共6个.而根据文章中构造的经验分布函数,13.018这个分界点对应的概率为0.9813,大于0.98.尽管0.98和0.9813相差不大,但是仍然有可能遗漏符合要求(即98%分位数及以上)的数据.在这里的表2可能存在这种疏漏,而实际上研究者侥幸逃脱这种疏漏,98%分位数及以上的数据恰好就只有6个.试想一下,如果样本量是3190的话,0.98和0.9813小小之差就会遗漏符合要求的4个数据(即在0.98和0.9813之间的数据),实在是不严谨,没有严格做到“挑出每一种重金属污染浓度最高的2%”.
    • 文章中的问题并不仅仅在于表2,表1和表1下面的经验分布函数也出现错误:首先看一下表1,表1的频率中0.0003和0.0006都应该是0.003和0.006,小数点后多了一个0.然后看一下表1下的“经验分布函数”,第4行开始的概率都随之而错,正确值应该是0.9812,0.9843,0.9906,0.9906,0.9969,0.9969,1.
    • 还有一个不小的问题.在表1下构造的经验分布函数中,每一段区间的(累计)概率其实上接近区间中最右端的样本数据的(教材)经验分布函数累计概率,这样的话在选取“最高的2%”的时候就会造成困扰:区间[10.116,13.018)对应的文中的经验分布函数为0.9812,大于0.98喔,选取“最高的2%”的时候是不是应该把10.116当作截点呢?而实际上这样做的话依照(已经修改过的)表1,选取到的是0.956分位数以上的数据,不符合要求.为什么会这样?我觉得就是文章构造的“经验分布函数”的先天缺陷了:构造文中的“经验分布函数”的过程损耗了有用信息——原本好好的数值变量为什么要分割成分类变量呢?——这样的构造不仅没有对研究有帮助,反而会造成困扰.

(3)文章方法的改进

  • 其实改进很简单,直接挑出每一种重金属污染浓度最高的2%的数据进行重点分析就好.根本没有必要构造文中的“经验分布函数”.

第三题：

阅读课本P29页实验1.8,参照程序给出waiting变量的经验分布函数置信区间(自己编程制图并分析),说明置信度与置信区间的关系.

library(MASS)
data(geyser)
waiting<-geyser[,1]
###开始画经验分布函数
waiting.sort<-sort(waiting)
waiting.rank<-rank(waiting.sort,ties.method="max")
waiting.cdf<-waiting.rank/length(waiting)
plot(waiting.sort,waiting.cdf)#描点
N<-length(waiting)
segments(waiting.sort[1:(N-1)],waiting.cdf[1:(N-1)],waiting.sort[2:N],waiting.cdf[1:(N-1)])#画线
###开始画经验分布函数的置信区间
alpha<-0.5
band<-sqrt(1/(2*N))*log(2/alpha)
lower.95<-waiting.cdf-band
upper.95<-waiting.cdf+band
lines(waiting.sort,lower.95,lty=2)
lines(waiting.sort,upper.95)

###对比不同alpha值下的置信区间宽度
alpha<-seq(0.01,0.1,by=0.001)
band<-sqrt(1/(2*N))*log(2/alpha)
plot(alpha,band)

可见置信区间的宽度随α值的增大而变小.

  • 附:用ecdf函数画经验密度函数(添加置信区间同理可得)
###画经验密度函数
Fn <- ecdf(waiting)
plot(Fn)
###画置信区间
lines(waiting.sort,lower.95,lty=2)
lines(waiting.sort,upper.95)