【非参数统计】2.1广义符号检验（R语言）

最新推荐文章于 2024-07-11 02:08:11 发布

一只小小夕

最新推荐文章于 2024-07-11 02:08:11 发布

阅读量7.6k

点赞数 15

分类专栏：非参数统计文章标签： r语言算法

本文链接：https://blog.csdn.net/weixin_50466307/article/details/115901005

版权

非参数统计专栏收录该内容

1 篇文章 4 订阅

订阅专栏

2.1 广义符号检验

1、符号检验结果（样本数据应支持备择假设）

法1：自动识别零假设（根据样本分位数大小，设为备择假设）

##自动识别零假设
sign.test=function(x,p,M0)   #x为数据，p为分位数，M0为待检验的的数
{s1=sum(x<M0);s2=sum(x>M0);n=s1+s2
p1=pbinom(s1,n,p);p2=1-pbinom(s1-1,n,p)
if (p1>p2) m1="H0: M>=M0"
else m1="H0: M<=M0"
p.value=min(p1,p2);p.value2=2*p.value
list(c("s+"=s2,"s-"=s1,"n'"=n),c("原假设"=m1,"单边p值"=p.value,"双边p值"=p.value2))
}

法2：自己设立零假设（根据样本分位数大小，设为备择假设）
H $_0$ : M $\leqq$ M $_0$
H $_0$ : M $\geqq$ M $_0$
H $_0$ : M $=$ M $_0$

#H0: M<=M0
sign.test1=function(x,p,M0) 
{s1=sum(x<M0)
s2=sum(x>M0)
n=s1+s2
p.value=pbinom(s1,n,p)
list(c("s+"=s2,"s-"=s1,"n'"=n),"p值"=p.value)
}

#H0: M>=M0
sign.test2=function(x,p,M0)
{s1=sum(x<M0)
s2=sum(x>M0)
n=s1+s2
p.value=1-pbinom(s1-1,n,p)
list(c("s+"=s2,"s-"=s1,"n'"=n),"p值"=p.value)
}

#H0: M=M0
sign.test3=function(x,p,M0)
{s1=sum(x<M0)
s2=sum(x>M0)
n=s1+s2
p.value1=1-pbinom(s1-1,n,p)
p.value2=pbinom(s1,n,p)
p.value=2*min(p.value1,p.value2)
list(c("s+"=s2,"s-"=s1,"n'"=n),"双边p值"=p.value)
}

- 例题：71个大城市的花费指数如下

74.5 74.3 73.9 71.7 71.2 67.7 66.7 66.2 65.4 65.3 65.3 65.3 64.6 63.5
62.7 60.8 58.2 55.5 55.3 55 54.9 52.7 51.8 49.9 48.2 47.6 46 45.8 45.2
41.9 38.8 37.7 37.5 36.5 36.4 32.7 32.7 32.2 29.1 27.8 27.8
请检验以下问题：
（1）有人说64是中位数（样本中位数为67.7，大于64）
（2）有人说64是下四分位数（样本下四分位数为50.85，小于64）

setwd("E:/非参数统计/data")

其次，导入数据（此处为txt数据，数据无标题行）

x1<-read.table("ExpensCities.txt",sep=" ")

法1：自动识别零假设

sign.test(x1,0.5,64)
sign.test(x1,0.25,64)

结果如下：

> sign.test(x1,0.5,64)
[[1]]
s- s+ n' 
28 43 71 

[[2]]
              原假设              单边p值              双边p值 
         "H0: M<=M0" "0.0479618157054121" "0.0959236314108241"

结论：s $^-$ =28，s $^+$ =43，n’=s $^-$ +s $^+$ =71，p 值=0.04796182<0.05

$\quad\quad\;$ 在 $\alpha$ =0.05时，拒绝原假设，认为中位数大于64

[[1]]
s- s+ n' 
28 43 71 

[[2]]
               原假设               单边p值               双边p值 
          "H0: M>=M0" "0.00515187949200158"  "0.0103037589840032"

结论：s $^-$ =28，s $^+$ =43，n’=s $^-$ +s $^+$ =71，p 值=0.005151879<0.01

$\quad\quad\;$ 在 $\alpha$ =0.01时，拒绝原假设，认为下四分位数小于64

法2：自己设立零假设

由于样本中位数为67.7，大于64，设 H $_0$ : M $\leqq$ M $_0$

sign.test1(x1,0.5,64)

结果如下：

> sign.test1(x1,0.5,64)
[[1]]
s- s+ n' 
28 43 71 

$p值
[1] 0.04796182

结论：s $^-$ =28，s $^+$ =43，n’=s $^-$ +s $^+$ =71，p 值=0.04796182<0.05

$\quad\quad\;$ 在 $\alpha$ =0.05时，拒绝原假设，认为中位数大于64

由于样本下四分位数为50.85，小于64，设 H $_0$ : M $\geqq$ M $_0$

sign.test2(x1,0.25,64)

结果如下：

> sign.test2(x1,0.25,64)
[[1]]
s- s+ n' 
28 43 71 

$p值
[1] 0.005151879

结论：s $^-$ =28，s $^+$ =43，n’=s $^-$ +s $^+$ =71，p 值=0.005151879<0.01

$\quad\quad\;$ 在 $\alpha$ =0.01时，拒绝原假设，认为下四分位数小于64

2、符号检验置信区间

qci=function(x,alpha=0.05,p=.25){         #x为数据,alpha为置信度,p为分位数
  x<-sort(x);n=length(x);a=alpha/2;r=qbinom(a,n,p);
  s=qbinom(1-a,n,p);CL=pbinom(s,n,p)-pbinom(r-1,n,p)
  if (r==0) lo<-NA else lo<-x[r]
  if (s==n) up<-NA else up<-x[s+1]
  list(c("置信下限"=lo,"置信上限"=up,
         "置信度"=1-alpha,"实际置信度"=CL),c("置信下限位置"=r,"置信上限位置"=s+1)) 
}

注：数据x必须先排序（sort()函数）

- 例题：同上，求花费指数中位数95%的置信区间

qci(sort(x1[,1]),0.05,0.5)

结果如下：

> qci(sort(x1[,1]),0.05,0.5)
[[1]]
  置信下限   置信上限     置信度 实际置信度 
62.7000000 77.7000000  0.9500000  0.9680728 

[[2]]
置信下限位置 置信上限位置 
          27           45

结论：花费指数中位数的置信区间[62.7,77.7]，实际置信度为96.8%

$\quad\;\;\;\;\;$ 置信位置[x $_($ $_2$ $_7$ $_)$ ,x $_($ $_4$ $_5$ $_)$ ]

2.2 cox-stuart趋势存在性检验

由于n为奇数偶数时，数对个数不一样。因此，以下函数分类表述：

#n是偶数
cox.test1=function(x)   #x为数据，p为分位数，M0为待检验的的数
{n=nrow(x)
D=x[1:(n/2),]-x[(n/2+1):n,]
s1=sum(D<0);s2=sum(D>0);n1=s1+s2;
p.value1=pbinom(min(s1,s2),n1,0.5)
p.value2=2*p.value1
if (s1>s2) m1="H0: 没有上升趋势"
else m1="H0: 没有下降趋势"
list(c("s-"=s1,"s+"=s2,"n'"=n1),c("原假设"=m1,"单边p值"=p.value1,"双边p值"=p.value2))
}
#n是奇数
cox.test2=function(x)   #x为数据，p为分位数，M0为待检验的的数
{n=nrow(x)
D=x[1:((n+1)/2-1),]-x[((n+1)/2+1):n,]
s1=sum(D<0);s2=sum(D>0);n1=s1+s2;
p.value1=pbinom(min(s1,s2),n1,0.5)
p.value2=2*p.value1
if (s1>s2) m1="H0: 没有上升趋势"
else m1="H0: 没有下降趋势"
list(c("s-"=s1,"s+"=s2,"n'"=n1),c("原假设"=m1,"单边p值"=p.value1,"双边p值"=p.value2))
}

- 例题：天津108个月的旅客数如下，问有无趋势？

54379 45461 55408 59712 60776 57635 63335 71296 70250 76866 75561 66427 61330 58186 67799 76360 86207 75509 83020 89614 75791 80835 72179 61520 66726 60629 68549 73310 80719 67759 70352 82825 70541 74631 68938 53318 62653 58578 63292 69535 73379 62859 72873 87260 67559 76647 70590 58935 58161 64057 63051 58807 63663 57367 70854 79949 66992 80140 62260 55942
58367 56673 61039 74958 85859 67263 87183 97575 79988 88501 68600 58442 68955 56835 67021 81547 85118 70145 95080 106186 86103 88548 70090 65550 69223 85138 89799 99513 98114 68172 97366 116820 95665 109881 87068 75362 88268 85183 87909 79976 27687 50178 100878 131788 116293 120770 104958 109603

由于n=108为偶数，因此选择cox.test1函数进行检验
而样本中s $^-$ =38，s $^+$ =16，负号较多，表明有增长的趋势
H $_0$ : M $\leqq$ M $_0$

x2<-read.csv("TJAir.csv",header=F,sep="")
cox.test1(x2)

结果如下：

> x2<-read.csv("TJAir.csv",header=F,sep="")
> cox.test1(x2)
[[1]]
s- s+ n' 
38 16 54 

[[2]]
               原假设               单边p值               双边p值 
   "H0: 没有上升趋势" "0.00191913294016305"  "0.0038382658803261"

结论：

附录：完整代码

#####################第2章 单样本位置检验############################
setwd("E:/武慧丽/大学/大三下/非参数统计/data")

##2.1 广义符号检验---------------------------------------------------

##< 1 符号检验结果
##自动识别零假设
sign.test=function(x,p,M0)   #x为数据，p为分位数，M0为待检验的的数
{s1=sum(x<M0);s2=sum(x>M0);n=s1+s2
p1=pbinom(s1,n,p);p2=1-pbinom(s1-1,n,p)
if (p1>p2) m1="H0: M>=M0"
else m1="H0: M<=M0"
p.value=min(p1,p2);p.value2=2*p.value
list(c("s+"=s2,"s-"=s1,"n'"=n),c("原假设"=m1,"单边p值"=p.value,"双边p值"=p.value2))
}

#H0: M<=M0
sign.test1=function(x,p,M0) 
{s1=sum(x<M0)
s2=sum(x>M0)
n=s1+s2
p.value=pbinom(s1,n,p)
list(c("s+"=s2,"s-"=s1,"n'"=n),"p值"=p.value)
}

#H0: M>=M0
sign.test2=function(x,p,M0)
{s1=sum(x<M0)
s2=sum(x>M0)
n=s1+s2
p.value=1-pbinom(s1-1,n,p)
list(c("s+"=s2,"s-"=s1,"n'"=n),"p值"=p.value)
}

#H0: M=M0
sign.test3=function(x,p,M0)
{s1=sum(x<M0)
s2=sum(x>M0)
n=s1+s2
p.value1=1-pbinom(s1-1,n,p)
p.value2=pbinom(s1,n,p)
p.value=2*min(p.value1,p.value2)
list(c("s+"=s2,"s-"=s1,"n'"=n),"双边p值"=p.value)
}

#书例2.1-P16
x1<-read.table("ExpensCities.txt",sep=" ")
sign.test(x1,0.5,64)
sign.test(x1,0.25,64)

sign.test1(x1,0.5,64)
sign.test2(x1,0.25,64)

## < 2 基于符号检验分位数的置信区间
qci=function(x,alpha=0.05,p=.25){         #x为数据,alpha为置信度,p为分位数
  x<-sort(x);n=length(x);a=alpha/2;r=qbinom(a,n,p);
  s=qbinom(1-a,n,p);CL=pbinom(s,n,p)-pbinom(r-1,n,p)
  if (r==0) lo<-NA else lo<-x[r]
  if (s==n) up<-NA else up<-x[s+1]
  list(c("置信下限"=lo,"置信上限"=up,
         "置信度"=1-alpha,"实际置信度"=CL),c("置信下限位置"=r,"置信上限位置"=s+1)) 
}

qci(sort(x1[,1]),0.05,0.5)