背景
假设我们有一个字母表A,B, C, D,然后我们查看一些数据并找到一个“单词”,DDDDDDDDCDDDDDD在我看来,找到这个随机的机会似乎很低,而发现BABDCABCDACDBACD似乎不太随机。
问题
我应该如何检查我遇到的字符串是否不是随机的?
我在 R 中尝试了一些东西,例如,对字母进行数字编码,然后将它们与排列进行比较。但是预先编码是相当麻烦的。可能有更直接的方法吗?
背景
假设我们有一个字母表A,B, C, D,然后我们查看一些数据并找到一个“单词”,DDDDDDDDCDDDDDD在我看来,找到这个随机的机会似乎很低,而发现BABDCABCDACDBACD似乎不太随机。
问题
我应该如何检查我遇到的字符串是否不是随机的?
我在 R 中尝试了一些东西,例如,对字母进行数字编码,然后将它们与排列进行比较。但是预先编码是相当麻烦的。可能有更直接的方法吗?
找到这个随机的机会对我来说似乎很低,而找到 BABDCABCDACDBACD 似乎不太随机。
为什么会这样?如果每个字母 A...D 的总比例等于 0.25,并且每个字母独立于另一个,那么这两个词的概率完全相同。如果字母的分布不同,那么生成两个单词的概率当然可能不同。
您可以尝试找到“低复杂度”的单词,例如一个字母比例特别高的单词(您可以使用其他响应中建议的香农信息,在生物序列分析中还有许多其他方法),但是不是对“随机性”的测试,因为没有进一步的假设或对您实际分析的内容的了解,“随机性”一词毫无意义。
你可以试试香农信息:
对于你的第一个词. 在第二个词中你有.
如果熵很高,你可以认为它比另一个熵较低的词更随机。
这里的其他答案集中在序列中不同字母的整体出现上,这可能是预期的“随机性”的一个方面。然而,另一个令人感兴趣的方面是序列中字母顺序的明显随机性。至少,我认为“随机性”需要字母向量的可交换性,可以使用“运行测试”进行测试。运行测试计算序列中“运行”的数量,并将运行总数与其在可交换性零假设下的零分布进行比较,对于具有相同字母的向量。构成“运行”的确切定义取决于特定的测试(例如,请参见此处的类似答案),但在这种情况下,对于名义类别,自然定义是将仅包含一个字母的任何连续序列计为单个“运行”。
例如,你的序列在我BABD-CABC-DACD-BACD看来表面上是非随机的(没有字母出现在它本身上,这对于这么长的序列来说可能是不寻常的)。 为了正式测试这一点,我们可以进行可交换性运行测试。在这个序列中,我们有字母(每个字母四个)并且有运行,每个由一个字母的单个实例组成。在可交换性假设下,可以将观察到的运行次数与其零分布进行比较。我们可以通过模拟来做到这一点,这会产生一个模拟的零分布和一个用于测试的 p 值。此字符序列的结果如下图所示。
对于这个序列,游程检验的 p 值(在可交换性的原假设下)是. 这在 10% 的显着性水平上是显着的,但在 5% 的显着性水平上并不显着。有一些证据表明存在不可交换的序列(即非随机顺序),但证据并不是特别有力。使用更长的观察字符串,运行测试将具有更大的能力来区分可更换的字符串和不可更换的字符串。(如您所见,我最初初步判断该字符串是非随机的可能是错误的——p 值实际上并没有我预期的那么低。)
最后,需要注意的是,这个测试只关注字符串中字母顺序的随机性——它将每种类型的字母数量作为固定输入。该测试将检测字符串中字母不可交换性意义上的非随机性,但不会测试不同字母整体概率意义上的“随机性”。如果后者也是“随机性”的指定含义的一部分,那么这个运行测试可以增加另一个测试,该测试查看字母的总计数,并将其与假设的零分布进行比较。
R 代码:上面的图和 p 值是使用以下R代码生成的:
#Define the character string vector (as factors)
x <- factor(c(2,1,2,4, 3,1,2,3, 4,1,3,4, 2,1,3,4),
labels = c('A', 'B', 'C', 'D'))
#Define a function to calculate the runs for an input vector
RUNS <- function(x) { n <- length(x);
R <- 1;
for (i in 2:n) { R <- R + (x[i] != x[i-1]) }
R; }
#Simulate the runs statistic for k permutations
k <- 10^5;
set.seed(12345);
RR <- rep(0, k);
for (i in 1:k) { x_perm <- sample(x, length(x), replace = FALSE);
RR[i] <- RUNS(x_perm); }
#Generate the frequency table for the simulated runs
FREQS <- as.data.frame(table(RR));
#Calculate the p-value of the runs test
R <- RUNS(x);
R_FREQ <- FREQS$Freq[match(R, FREQS$RR)];
p <- sum(FREQS$Freq*(FREQS$Freq <= R_FREQ))/k;
#Plot estimated distribution of runs with test
library(ggplot2);
ggplot(data = FREQS, aes(x = RR, y = Freq/k, fill = (Freq <= R_FREQ))) +
geom_bar(stat = 'identity') +
geom_vline(xintercept = match(R, FREQS$RR)) +
scale_fill_manual(values = c('Grey', 'Red')) +
theme(legend.position = 'none',
plot.title = element_text(hjust = 0.5, face = 'bold'),
plot.subtitle = element_text(hjust = 0.5),
axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0))) +
labs(title = 'Runs Test - Plot of Distribution of Runs under Exchangeability',
subtitle = paste0('(Observed runs is black line, p-value = ', p, ')'),
x = 'Runs', y = 'Estimated Probability');
我已经用破折号打破了序列,只是为了更容易阅读;破折号对分析没有意义。