全部版块 我的主页
论坛 数据科学与人工智能 数据分析与数据科学 R语言论坛
17048 7
2014-04-04
对一个序列做单样本K-S检验 检验其是否符合[0,1]均匀分布(此序列是由一组时间序列进行计算累积分布概率得出的序列)
检验使用的语句如下:ks.test(a,"punif",0,1)
  检验结果是
1.jpg
p值特别小,不符合[0,1]均匀分布

但是对这个序列做[0,1]均匀分布的Q-Q图,却几乎完全拟合
2.jpg

请问:检验结果和图形拟合不同,这是为什么啊? 到底是哪里出错了?该怎么做才能将一个给定的时间序列转化成[0,1]均匀分布?



二维码

扫码加我 拉你入群

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

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

全部回复
2014-4-4 20:57:25
可能序列有一定的自相关性
二维码

扫码加我 拉你入群

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

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

2014-10-9 21:50:10
R中的ks.test检验确实有点不靠谱。也许是个人水平不够,该函数用起来非常不顺手。我主要用wilcox.test来检验均匀分布。举个例子:先随机产生500个数,每个数的取值范围是:大于等于1且小于等于10,这500个数都为整数。现在要求检验这500个数是否符合均匀分布。假设这500个数存于向量x中,那么,产生一组与其长度差不多的向量y,其中也同为整数,且依据均匀分布产生1到10之间的数字。之后用wilcox.test(x,y)检测两组数据的均值是否存在显著不同(p值小于临界值则认为不同),由于该函数不事先假设样本的分布,所以更贴近实际情况。但是,光做一次比较不能说明问题:即使此次检验p值很大,也不能说x就服从均匀分布。均匀分布有个特点:

比如均匀分布产生的向量z=(2,1,3,4,5,9,8,6,7,10, ...) ---->将所有数字全部加1,超过10的再用1替代,形成“环形位移”---->形成新的向量z' = (3,2,4,5,6,10,9,7,8,1, ...)   :  z与z' 的均值应该都差不多!当然,在这个例子中两者是相等的,实际例子中会有些许差异,但相差不会太大,原因是——总共10个档位(取值范围是1-10),无论你怎么交换每个档位代表的数值大小,每个档位上聚集的数字个数彼此都差不多。比如500个数字,分到10个档位,每个档平均有500/10 = 50个左右的数字。因此,就算是在z' 的基础上再进行“环形位移”,产生z'',并在此基础上再位移,产生z''',它们彼此的均值水平都相差无几!

如果向量z并不服从均匀分布呢?比如z = (3,3,3,3,3,3,3,3),那么位移后为z' = (4,4,4,4,4,4,4,4),均值差异明显,更不要说拿几轮位移之后的z'''和z去比较了。

因此,检验均匀分布,除了运行wilcox.test(x,y),还要运行wilcox.test(x',y), wilcox.test(x'',y), wilcox.test(x''',y), ...位移要进行9次(档位总共只有10档,第10次位移就又变回向量x了,所以没有必要再做了),检验也要多进行9次,再加上第一次检验,如果10次中有比如9次能通过(理想情况是10次全过),则说明向量x中的数据很大可能性是来自均匀分布!

程序如下:

# mov为环形位移程序(是unif.test的子函数),参数列表如下:
# x : 需要位移的数据向量
# maxValue : x中数据若加1后超过该数,就会调成最小值(例子中maxValue为10,最小值系统会自动计算,此处为1)
# n : 档位数,由于每个整数随机变量取值范围是1-10,故档位总共有10个
# 函数会返回一个与x等长的向量,它是环形位移1位之后的结果

mov <- function(x, maxValue, n) {

  x = x + 1
  x[which(x > maxValue)] = x[which(x > maxValue)] - n

  return(x)

}

# 均匀分布检验程序unif.test,参数列表如下:
# x : 需进行均匀分布检验的数据向量
# n : 档位数,调用函数时若第二位不写,则默认为10
# alpha : 判断p值是否过小的临界值,调用函数时若第三位不写,则默认为0.1,p值小于该数,则视为x均值与y显著不同
# 函数返回一个小数,该小数是用n次测试中成功的次数(即x或其位移向量与y均值大体相同)除以总测试次数n所得到的比例
# 按上例,该数应大于等于0.9,即10次中至少要有9次成功才可以,否则就太不像话了,不能承认x是服从均匀分布的
# 函数用法:unif.test(x) ---> 只给出待检验向量,后面空着,则n默认为10,alpha默认为0.1,当然,也可赋其他值,比如:
#                   unif.test(x, 100, 0.05),但在题目所给情形下这样定参数是没道理的

unif.test <- function(x, n = 10, alpha = 0.1) {

  max.x = max(x)
  rep.y = round(length(x) / n)
  if (rep.y < 2) rep.y = 2
  y = rep((max.x - n + 1):max.x, rep.y)   # -------> 产生标准样本y

  result = 0

  for (i in 1:n) {
  
    x = mov(x, max.x, n)
    wt = wilcox.test(x, y)   # -----> 默认双侧检验

    if (wt$p.value > alpha) result = result + 1
  
  }

  return(result/n)

}


测试:

> x = ceiling(runif(500)*10)  # --------> 产生500个整数随机变量,每个的取值范围是1-10,这里用均匀分布产生随机数
> unif.test(x)
[1] 1

结果是合格的,数据符合均匀分布

如果产生整数随机数时使用正态分布,那么检验结果如何?

> z = seq(from = -3, to = 3, by = 0.6)
> p = diff(pnorm(z))
> p[5] = p[5] + pnorm(-3)
> p[6] = p[6] + pnorm(-3)
> x = sample(1:10, 500, replace = TRUE, prob = p)   # -------> 其中的p表示1到10出现的概率,它是正态分布对应的概率
> unif.test(x)
[1] 0.2

结果明显不合格,样本不符合均匀分布

二维码

扫码加我 拉你入群

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

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

2014-10-10 20:43:37
OK,接上贴。以上均匀分布检验仅处理了一组500个数据的情况,如果是小样本呢?比如样本x仅含7个整数,检验会有效么?

况且就算是针对500个数据,以上实际运行也仅进行了一次,unif.test(x)的判断基本正确。但是这种检验的稳定性怎样?如果我重复做1000次这样的检验,都能判断正确么?

为此,我们考察样本容量从7变到100时每种情况的检验,针对其中的每个特定样本容量,各做1000次检验(即产生1000组样本容量为n个整数的样本,对它们进行unif.test检验),于是针对每个特定样本容量,都会得到1000个成功比例数,采集这1000个结果的:期望值、中位数、结果等于1的次数、结果大于等于0.9的次数、结果大于等于0.8的次数、结果大于等于0.7的次数、结果大于等于0.6的次数、结果大于等于0.5的次数、1000个结果中的最小值。


(1)如果我们用均匀分布产生样本,并对其作均匀分布检验unif.test,其结果如下:

i.range = 7:100
resultMatrix = matrix(nrow = length(i.range), ncol = 9)
colnames(resultMatrix) = c("mean", "median", "EQ1", "GE0.9", "GE0.8", "GE0.7", "GE0.6", "GE0.5", "min")
rownames(resultMatrix) = paste("X", i.range, sep = "")
row.index = 1
m.row = 1000

for (i in i.range) {

  m = matrix(nrow = m.row, ncol = i)
  for (j in 1:m.row) {
    m[j, ] = ceiling(runif(i)*10)
    m[j, m[j, ] == 0] = 1
  }
  result = as.vector(apply(m, 1, unif.test))

  resultMatrix[row.index, 1] = mean(result)
  resultMatrix[row.index, 2] = median(result)
  for (j in 0:5) resultMatrix[row.index, j + 3] = sum(result >= (1-j/10))
  resultMatrix[row.index, 9] = min(result)
  row.index = row.index + 1
  
}

> resultMatrix
       mean median EQ1 GE0.9 GE0.8 GE0.7 GE0.6 GE0.5 min
X7   0.9556      1 699   899   967   991  1000  1000 0.6
X8   0.9632      1 754   912   975   991  1000  1000 0.6
X9   0.9650      1 751   927   978   996   998  1000 0.5
X10  0.9605      1 743   911   971   985   995  1000 0.5
X11  0.9662      1 767   933   974   991   998   999 0.4
X12  0.9710      1 806   934   977   993  1000  1000 0.6
X13  0.9659      1 763   924   980   994   998  1000 0.5
X14  0.9716      1 802   943   983   992   996  1000 0.5
X15  0.9737      1 814   948   982   994   999  1000 0.5
X16  0.9717      1 795   945   984   994   999  1000 0.5
X17  0.9804      1 862   960   985   997  1000  1000 0.6
X18  0.9816      1 860   964   992  1000  1000  1000 0.7
X19  0.9769      1 849   944   984   994   998  1000 0.5
X20  0.9838      1 875   972   992   999  1000  1000 0.6
... ... ... ... ... ... ... ...
X90  0.9785      1 856   955   982   994   998  1000 0.5
X91  0.9804      1 862   964   988   995   997   998 0.4
X92  0.9820      1 866   967   989   998  1000  1000 0.6
X93  0.9815      1 852   971   993   999  1000  1000 0.6
X94  0.9844      1 875   978   993   998  1000  1000 0.6
X95  0.9766      1 833   953   983   997  1000  1000 0.6
X96  0.9778      1 832   962   988   997   999  1000 0.5
X97  0.9811      1 855   973   988   997   999   999 0.4
X98  0.9834      1 874   972   991   997  1000  1000 0.6
X99  0.9835      1 871   971   994   999  1000  1000 0.6
X100 0.9834      1 868   973   995   998  1000  1000 0.6

基本上对每种样本容量,unif.test检验结果在0.9或以上的概率均超90%,即1000次检验,有超900次其结果大于等于0.9。在样本容量大于等于17的情况下,其概率基本保持在95%以上。

至少针对1-10整数随机均匀分布的情况,unif.test比例底限确定在0.9,均匀分布被成功检测出的概率超90%,大多数情况超95%。


(2)如果我们用正态分布产生样本(注意,是正态分布,且同样是1-10范围内的整数随机样本),并对其作均匀分布检验unif.test,其结果如下:

z = seq(from = -3, to = 3, by = 0.6)
p = diff(pnorm(z))
p[5] = p[5] + pnorm(-3)
p[6] = p[6] + pnorm(-3)

i.range = 7:100
resultMatrix = matrix(nrow = length(i.range), ncol = 9)
colnames(resultMatrix) = c("mean", "median", "EQ1", "GE0.9", "GE0.8", "GE0.7", "GE0.6", "GE0.5", "min")
rownames(resultMatrix) = paste("X", i.range, sep = "")
row.index = 1
m.row = 1000

for (i in i.range) {

  m = matrix(nrow = m.row, ncol = i)
  for (j in 1:m.row) m[j, ] = sample(1:10, i, replace = TRUE, prob = p)

  result = as.vector(apply(m, 1, unif.test))

  resultMatrix[row.index, 1] = mean(result)
  resultMatrix[row.index, 2] = median(result)
  for (j in 0:5) resultMatrix[row.index, j + 3] = sum(result >= (1-j/10))
  resultMatrix[row.index, 9] = min(result)
  row.index = row.index + 1
  
}

> resultMatrix
       mean median EQ1 GE0.9 GE0.8 GE0.7 GE0.6 GE0.5 min
X7   0.7929    0.8 175   387   625   784   962   996 0.4
X8   0.7847    0.8 167   378   589   762   959   992 0.4
X9   0.7744    0.8 128   340   580   780   925   991 0.4
X10  0.7311    0.7 130   302   473   649   833   930 0.3
X11  0.7325    0.7 115   257   456   672   876   950 0.3
X12  0.7262    0.7  94   233   446   658   875   959 0.3
X13  0.6979    0.7  71   194   367   597   830   924 0.3
X14  0.6773    0.7  46   159   324   537   784   927 0.3
X15  0.6838    0.7  47   154   343   558   796   943 0.3
X16  0.6598    0.7  28   127   286   519   755   891 0.3
X17  0.6711    0.7  40   136   320   530   765   925 0.3
X18  0.6587    0.6  34   130   278   495   751   905 0.3
X19  0.6369    0.6  34   110   254   417   674   889 0.3
X20  0.6356    0.6  25   112   249   424   669   884 0.3
X21  0.6235    0.6  23    97   212   396   659   859 0.3
X22  0.6102    0.6  26    77   178   373   619   839 0.2
X23  0.6053    0.6  11    62   164   355   617   851 0.3
X24  0.6024    0.6  12    57   164   352   612   843 0.2
X25  0.5935    0.6  12    56   163   322   566   831 0.3
X26  0.5181    0.5   2    22    60   149   355   660 0.2
X27  0.5098    0.5   7    20    46   138   353   629 0.2
X28  0.5187    0.5   2    16    57   158   381   649 0.2
X29  0.4998    0.5   2    14    37   114   307   590 0.2
X30  0.4847    0.5   7    17    48   111   275   530 0.2
X31  0.4912    0.5   1     8    28   101   296   564 0.2
X32  0.4821    0.5   1     9    36   105   270   520 0.2
X33  0.4821    0.5   1     7    23    83   277   527 0.2
X34  0.4745    0.5   2     7    22    80   239   528 0.2
X35  0.4281    0.4   0     2     9    39   125   344 0.2
X36  0.4273    0.4   0     1    10    37   128   327 0.2
X37  0.4155    0.4   1     2     3    26   107   306 0.2
X38  0.4159    0.4   0     2     5    20   106   304 0.2
X39  0.4177    0.4   0     0     6    27   106   304 0.2
X40  0.4171    0.4   0     1     3    26   115   308 0.2
X41  0.4024    0.4   0     0     3    14    87   253 0.2
X42  0.4131    0.4   0     1     4    27    92   299 0.2
X43  0.4013    0.4   0     0     1    10    83   263 0.2
X44  0.3959    0.4   0     0     1    13    68   234 0.2
X45  0.3935    0.4   0     2     4    12    62   220 0.2
X46  0.3607    0.4   0     0     0     1    23   144 0.2
X47  0.3556    0.3   0     0     2     5    26   127 0.1
X48  0.3590    0.4   0     0     0     2    35   130 0.1
X49  0.3534    0.3   0     0     0     5    26   105 0.1
X50  0.3472    0.3   0     0     0     1    23   104 0.1
X51  0.3425    0.3   0     0     0     1    12    88 0.1
X52  0.3432    0.3   0     0     0     2    13    87 0.1
X53  0.3457    0.3   0     0     0     0    20   112 0.2
X54  0.3400    0.3   0     0     1     2    12    78 0.2
X55  0.3206    0.3   0     0     0     0    11    57 0.1
X56  0.3128    0.3   0     0     0     0     5    40 0.1
X57  0.3103    0.3   0     0     0     0     9    52 0.1
X58  0.3050    0.3   0     0     0     0     2    40 0.1
... ... ... ... ... ... ... ...
X95  0.2151    0.2   0     0     0     0     0     0 0.1
X96  0.2137    0.2   0     0     0     0     0     0 0.1
X97  0.2125    0.2   0     0     0     0     0     0 0.1
X98  0.2143    0.2   0     0     0     0     0     1 0.1
X99  0.2115    0.2   0     0     0     0     0     1 0.1
X100 0.2093    0.2   0     0     0     0     0     0 0.1

当样本容量较小时,即使unif.test结果的底限设定在0.9,其误诊概率还是很高:样本容量小于20的情况下,将例子中特定正态分布样本检验成均匀分布的误诊概率超过10%,最高可达38.7%,即1000次检验中,至少有超过100次,无法拒绝样本来自均匀分布(当然,这是不对的,因为我们是用正态分布的概率产生随机数的)。不过当样本容量大于35时,在底限为0.9的水平上,基本都可正确检验出样本并非来自均匀分布。其误诊概率几乎为0%,1000次检验中,最多也仅有几次误诊。


(3)作为对比,我们看一下ks.test的检验效果:首先,该函数不宜检验整型随机数是否服从均匀分布,效果不好(至少在整数取值为1-10这种“狭窄范围”的情况下,不宜用ks.test检验均匀分布)。其次,就算是检验连续随机变量的情况,其效果也不太理想。我们用均匀分布产生随机数,看看ks.test是否能检验出其服从均匀分布:

i.range = 7:100
resultVector = vector(length = length(i.range))
names(resultVector) = paste("X", i.range, sep = "")
v.index = 1

for (i in i.range) {

  result = 0

  for (j in 1:1000) {
  
    x = runif(i)*10  # ----> 此处妥协一下,不限定随机变量必须为整数,但范围与前例相近,取0到10。放弃整数是因为ks.test不宜做整数的均匀分布检验(至少在整数取值为1-10这种“狭窄范围”的情况下,不宜用ks.test检验均匀分布)
    ks = ks.test(x, "punif", min = 0, max = 10)
    if (ks$p.value > 0.1) result = result + 1
  
  }

  resultVector[v.index] = result
  v.index = v.index + 1

}

> resultVector
  X7   X8   X9  X10  X11  X12  X13  X14  X15  X16  X17  X18  X19  X20  X21
907  915  915  901  915  908  895  911  900  900  895  894  891  915  909
X22  X23  X24  X25  X26  X27  X28  X29  X30  X31  X32  X33  X34  X35  X36
906  899  902  913  895  904  901  901  905  902  899  898  905  896  891
X37  X38  X39  X40  X41  X42  X43  X44  X45  X46  X47  X48  X49  X50  X51
885  892  895  915  911  890  912  881  915  901  897  907  904  916  901
X52  X53  X54  X55  X56  X57  X58  X59  X60  X61  X62  X63  X64  X65  X66
893  888  913  919  902  911  894  904  896  900  901  914  893  911  907
X67  X68  X69  X70  X71  X72  X73  X74  X75  X76  X77  X78  X79  X80  X81
898  888  906  911  895  909  898  879  898  901  913  893  886  895  897
X82  X83  X84  X85  X86  X87  X88  X89  X90  X91  X92  X93  X94  X95  X96
913  914  898  897  910  893  888  910  901  900  888  900  917  897  886
X97  X98  X99 X100
910  908  891  905

样本容量仍从7变到100,每个特定容量水平,产生1000个样本,并记录ks.test成功检测出其服从均匀分布的次数(以p值大于0.1为标准)。结果如上所示。

首先说明,ks.test检测的稳定性不错,其结果不会依样本容量的扩大而产生明显变化。但是其检测成功的概率也就基本维持在90%左右,其整体表现不如unif.test。unif.test在样本容量充足的情况下,检验成功率保持在95%以上。
二维码

扫码加我 拉你入群

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

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

2015-3-28 13:24:53
我也是这个问题,你是怎么处理的?
二维码

扫码加我 拉你入群

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

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

2018-3-10 01:35:46
想问一下后来是怎么处理的哈?
二维码

扫码加我 拉你入群

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

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

点击查看更多内容…
相关推荐
栏目导航
热门文章
推荐文章

说点什么

分享

扫码加好友,拉您进群
各岗位、行业、专业交流群