##自动识别零假设
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)
}
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")
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:自己设立零假设
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
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
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()函数)
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 _2 2 7 _7 7 ) _) ),x ( _( ( 4 _4 4 5 _5 5 ) _) )]
由于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))
}
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
_0
0: M
≦
\leqq
≦ M
0
_0
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)
链接:https://pan.baidu.com/s/1Dh_odHA5LIBNYqElxkEGJg
提取码:lbjm