Brown-Mood检验 R程序

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

###Brown-Mood中位数检验（精确检验，正态近似，连续性修正后的正态近似）
BM.test=function(x,y,alt) #alt:备择假设形式
  {
  xy=c(x,y)
  md.xy=median(xy)
  t=sum(xy>md.xy)
  lx=length(x[x!=md.xy])
  ly=length(y[y!=md.xy])
  lxy=lx+ly
  A=sum(x>md.xy)#检验统计量A
  z=(A-lx*t)/(lx+ly)/(lx*ly*t*(lx+ly-t)/(lx+ly)^3)^0.5#正态近似时的标准化统计量
  if(A>(min(lx,t)/2)){
    z1=(A+0.5-lx*t)/(lx+ly)/(lx*ly*t*(lx+ly-t)/(lx+ly)^3)^0.5#连续性修正后正态近似时的标准化统计量
  }
  else{z1=(A-0.5-lx*t)/(lx+ly)/(lx*ly*t*(lx+ly-t)/(lx+ly)^3)^0.5}
  if(alt=="greater"){
    pv1=1-phyper(A,lx,ly,t)#精确p值
    pv2=1-pnorm(z)#正态近似p值
    pv3=1-pnorm(z1)#连续性修正后正态近似p值
  }
  if(alt=="less"){
    pv1=phyper(A,lx,ly,t)
    pv2=pnorm(z)
    pv3=pnorm(z1)
  }
  if(alt=="two.sided"){
    pv1=2*min(1-phyper(A,lx,ly,t),phyper(A,lx,ly,t))
    pv2=2*min(1-pnorm(z),pnorm(z))
    pv3=2*min(1-pnorm(z1),pnorm(z1))
  }
  conting.table=matrix(c(A,lx-A,lx,t-A,ly-(t-A),ly,t,lxy-t,lxy),3,3)#计数表
  col.name=c("X","Y","X+Y")
  row.name=c(">MXY","<MXY","TOTAL")
  dimnames(conting.table)=list(row.name,col.name)
  list(contingency.table=conting.table,p.value=pv1,pvnorm=pv2,pvnr=pv3)
}

###练习4.1
#Brown-Mood中位数检验
a=c(10,8,12,16,5,9,7,11,6)
b=c(12,15,20,18,13,14,9,16)
BM.test(a,b,alt="less")

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：Brown Mood Own R程序 function 程序 检验

运行错误，说没有BM。test这个函数。是不是需要装软件包啊？

额问问 发表于 2016-7-1 18:02
运行错误，说没有BM。test这个函数。是不是需要装软件包啊？

BM.test 这个函数是前面那些代码，这个函数是自己编写的。

# Brown-Mood检验##################################################
BMq.test = function(x,y,q,alt)
{
xy = c(x,y)
quantile.xy = quantile(xy,q)
t = sum(xy>quantile.xy)
lx = length(x[x!=quantile.xy])
ly = length(y[y!=quantile.xy])
lxy = lx + ly
A = sum(x>quantile.xy)        #检验统计量A
z = (A-lx*t/(lx+ly))/(lx*ly*t*(lx+ly-t)/(lx+ly)^3)^0.5
if(A>(lx*t/(lx+ly))){
        z1 = (A+0.5-lx*t/(lx+ly))/(lx*ly*t*(lx+ly-t)/(lx+ly)^3)^0.5
                                        #正态近似时的标准化统计量
}

else
{z1 = (A-0.5-lx*t/(lx+ly))/(lx*ly*t*(lx+ly-t)/(lx+ly)^3)^0.5}
if(alt == 'greater')
{
        pv1 = 1-phyper(A,lx,ly,t)
        pv2 = 1-pnorm(z)
        pv3 = 1-pnorm(z1)
}
if(alt == 'less')
{
        pv1 = phyper(A,lx,ly,t)
        pv2 = pnorm(z)
        pv3 = pnorm(z1)
}
if(alt == 'two.sided')
{
        pv1 = 2*min(1-phyper(A,lx,ly,t),phyper(A,lx,ly,t))
        pv2 = 2*min(1-pnorm(z),pnorm(z))
        pv3 = 2*min(1-pnorm(z1),pnorm(z1))
}
conting.table = matrix(c(A,lx-A,lx,t-A,ly-(t-A),ly,t,lxy-t,lxy),3,3)
col.name = c('X','Y','X+Y')
row.name = c('>MQXY','<MQXY','TOTAL')
dimnames(conting.table) = list(row.name,col.name)
list(contingency.table = conting.table,p.value = pv1,pvnorm = pv2,pvnr = pv3)
}
x = c(10,8,12,16,5,9,7,11,6)
y = c(12,15,20,18,13,14,9,16)
xy = c(x,y)
median(xy)
BMq.test(x,y,0.25,'two.sided')