@NickCox 提出了一个有趣的策略 (+1)。然而，由于@whuber指出的担忧，我可能认为它本质上更具探索性。

让我建议另一种策略：您可以拟合高斯有限混合模型。请注意，这是非常强烈的假设，即您的数据来自一个或多个真实的法线。正如@whuber 和@NickCox 在评论中指出的那样，如果没有对这些数据进行实质性解释（由完善的理论支持）来支持这一假设，该策略也应被视为探索性的。

首先，让我们按照@Glen_b 的建议，使用两倍的 bin 来查看您的数据：

我们仍然看到两种模式；如果有的话，他们在这里更清楚地出现了。（另请注意，内核密度线应该相同，但由于箱数较多而显得更加分散。）

现在让我们拟合一个高斯有限混合模型。在R中，您可以使用Mclust包来执行此操作：

library(mclust)
x.gmm = Mclust(x)
summary(x.gmm)
# ----------------------------------------------------
# Gaussian finite mixture model fitted by EM algorithm 
# ----------------------------------------------------
#   
# Mclust V (univariate, unequal variance) model with 2 components:
#   
#   log.likelihood   n df       BIC       ICL
#        -1200.874 120  5 -2425.686 -2442.719
# 
# Clustering table:
#  1  2 
# 68 52 

两个正常组件优化 BIC。为了比较，我们可以强制进行单分量拟合并执行似然比检验：

x.gmm.1 = Mclust(x, G=1)
logLik(x.gmm.1)
# 'log Lik.' -1226.241 (df=2)
logLik(x.gmm)-logLik(x.gmm.1)
# 'log Lik.' 25.36657 (df=5)
1-pchisq(25.36657, df=3)  # [1] 1.294187e-05

这表明，如果数据来自单个真正的正态分布，那么您极不可能找到与您的数据一样远离单峰的数据。

有些人在这里使用参数测试感到不舒服（尽管如果假设成立，我不知道有任何问题）。一种非常广泛适用的技术是使用参数引导交叉拟合方法（我在此处描述了该算法）。我们可以尝试将其应用于这些数据：

x.gmm$parameters
# $mean
# 12346.98 23322.06 
# $variance$sigmasq
# [1]  4514863 24582180
x.gmm.1$parameters
# $mean
# [1] 17520.91
# $variance$sigmasq
# [1] 43989870

set.seed(7809)
B = 10000;    x2.d = vector(length=B);    x1.d = vector(length=B)
for(i in 1:B){
  x2      = c(rnorm(68, mean=12346.98, sd=sqrt( 4514863)), 
              rnorm(52, mean=23322.06, sd=sqrt(24582180)) )
  x1      = rnorm( 120, mean=17520.91, sd=sqrt(43989870))
  x2.d[i] = Mclust(x2, G=2)$loglik - Mclust(x2, G=1)$loglik
  x1.d[i] = Mclust(x1, G=2)$loglik - Mclust(x1, G=1)$loglik
}
x2.d = sort(x2.d);  x1.d = sort(x1.d)
summary(x1.d)
#     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
# -0.29070 -0.02124  0.41460  0.88760  1.36700 14.01000 
summary(x2.d)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#  9.006  23.770  27.500  27.760  31.350  53.500 

抽样分布的汇总统计数据和核密度图显示了几个有趣的特征。单分量模型的对数似然很少大于两分量拟合的对数似然，即使真实数据生成过程只有一个分量，当它更大时，数量是微不足道的。比较具有不同拟合数据能力的模型的想法是 PBCM 背后的动机之一。两个采样分布几乎没有重叠；只有 0.35%x2.d小于最大值x1.d价值。如果您在对数似然差 > 9.7 时选择了双组件模型，那么您将错误地选择单组件模型 0.01% 和 0.02% 的双组件模型。这些是高度可区分的。另一方面，如果您选择使用一个组件模型作为原假设，那么您观察到的结果足够小，不会出现在 10,000 次迭代的经验抽样分布中。我们可以使用 3 规则（参见此处）为 p 值设置上限，即我们估计您的 p 值小于 0.0003。也就是说，这是非常重要的。

这就提出了一个问题，为什么这些结果与你的浸渍测试有很大的不同。（要回答您的明确问题，您的浸入测试没有提供任何证据表明存在两种真实模式。）老实说，我不知道浸入测试，因此很难说；它可能动力不足。但是，我认为可能的答案是这种方法假设您的数据是由真正的正常[s]生成的。对您的数据进行 Shapiro-Wilk 检验非常重要（$p < .000001$)，并且对于数据的最优 Box-Cox 变换（平方根倒数；$p < .001$）。但是，数据从来都不是真正正常的（参见这句著名的引言），并且如果存在底层组件，也不能保证它们完全正常。如果您发现您的数据可能来自正偏态分布而不是正态分布是合理的，那么这种双峰性水平很可能在典型的变化范围内，这就是我怀疑浸入测试所说的。

跟进@Nick 的回答和评论中的想法，您可以看到带宽需要多宽才能使辅助模式变平：

将此核密度估计作为最接近的零点（最接近数据但仍与零假设一致，即它是来自单峰总体的样本）并从中进行模拟。在模拟样本中，次要模式通常看起来并不那么明显，您无需尽可能扩大带宽以使其平坦。

将这种方法形式化导致 Silverman (1981)，“使用核密度估计来研究模态”，JRSS B， 43，1 中给出的测试。Schwaiger & Holzmann 的silvermantest包实现了这个测试，以及 Hall & York 描述的校准程序（ 2001），“关于 Silverman 多模态检验的校准”，Statistica Sinica，11，第 515 页，针对渐近保守主义进行了调整。使用单峰的零假设对您的数据执行检验会导致 p 值为 0.08（未经校准）和 0.02（经校准）。我对浸入测试不够熟悉，无法猜测为什么它可能会有所不同。

代码：

      # kernel density estimate for x using Sheather-Jones 
          # method to estimate b/w:
    density(x, kernel="gaussian", bw="SJ") -> dens.SJ
      # tweak b/w until mode just disappears:
    density(x, kernel="gaussian", bw=3160) -> prox.null
      # fill matrix with simulated samples from the proximal 
      # null:
    x.sim <- matrix(NA, nrow=length(x), ncol=10)
    for (i in 1:10){
      x.sim[ ,i] <- rnorm(length(x), sample(x, size=length(x), 
                          replace=TRUE), prox.null$bw)
    }
      # perform Silverman test without Hall-York calibration:
    require(silvermantest)
    silverman.test(x, k=1, M=10000, adjust=F)
      # perform Silverman test with Hall-York calibration:
    silverman.test(x, k=1, M=10000, adjust=T)

需要担心的事情包括：

1. 数据集的大小。它不小，也不大。

2. 您看到的内容对直方图原点和 bin 宽度的依赖性。只有一个选择显而易见，您（和我们）对敏感性一无所知。

3. 您所看到的内容对内核类型和宽度的依赖性以及在密度估计中为您做出的任何其他选择。只有一个选择显而易见，您（和我们）对敏感性一无所知。

在其他地方，我曾尝试性地建议，模式的可信度由实质性解释和在相同大小的其他数据集中辨别相同模式的能力来支持（但未建立）。（越大越好……）

我们不能在这里评论其中任何一个。关于可重复性的一个小技巧是比较你得到的相同大小的引导样本。以下是使用 Stata 进行令牌实验的结果，但您看到的内容被任意限制为 Stata 的默认设置，这些默认设置本身被记录为从空气中提取的。我得到了原始数据和 24 个自举样本的密度估计值。

我认为经验丰富的分析师会从您的图表中猜测出任何迹象（不多也不少）。左手模式是高度可重复的，右手模式显然更脆弱。

请注意，这是不可避免的：由于靠近右手模式的数据较少，它不会总是重新出现在引导样本中。但这也是重点。

请注意，上面的第 3 点保持不变。但结果介于单峰和双峰之间。

对于那些感兴趣的人，这是代码：

clear 
set scheme s1color 
set seed 2803 

mat data = (10346, 13698, 13894, 19854, 28066, 26620, 27066, 16658, 9221, 13578, 11483, 10390, 11126, 13487, 15851, 16116, 24102, 30892, 25081, 14067, 10433, 15591, 8639, 10345, 10639, 15796, 14507, 21289, 25444, 26149, 23612, 19671, 12447, 13535, 10667, 11255, 8442, 11546, 15958, 21058, 28088, 23827, 30707, 19653, 12791, 13463, 11465, 12326, 12277, 12769, 18341, 19140, 24590, 28277, 22694, 15489, 11070, 11002, 11579, 9834, 9364, 15128, 15147, 18499, 25134, 32116, 24475, 21952, 10272, 15404, 13079, 10633, 10761, 13714, 16073, 23335, 29822, 26800, 31489, 19780, 12238, 15318, 9646, 11786, 10906, 13056, 17599, 22524, 25057, 28809, 27880, 19912, 12319, 18240, 11934, 10290, 11304, 16092, 15911, 24671, 31081, 27716, 25388, 22665, 10603, 14409, 10736, 9651, 12533, 17546, 16863, 23598, 25867, 31774, 24216, 20448, 12548, 15129, 11687, 11581)
set obs `=colsof(data)' 
gen data = data[1,_n] 

gen index = . 

quietly forval j = 1/24 { 
    replace index = ceil(120 * runiform()) 
    gen data`j' = data[index]
    kdensity data`j' , nograph at(data) gen(xx`j' d`j') 
} 

kdensity data, nograph at(data) gen(xx d) 

local xstuff xtitle(data/1000) xla(10000 "10" 20000 "20" 30000 "30") sort 
local ystuff ysc(r(0 .0001)) yla(none) `ystuff'   

local i = 1 
local colour "orange" 
foreach v of var d d? d?? { 
    line `v' data, lc(`colour') `xstuff'  `ystuff' name(g`i', replace) 
    local colour "gs8" 
    local G `G' g`i' 
    local ++i 
} 

graph combine `G' 

LP Nonparametric Mode Identification（算法名称LPMode，论文参考如下）

MaxEnt模式 [图中的红色三角形]：12783.36 和 24654.28。

L2模式 [图中的绿色三角形]：13054.70 和 24111.61。

有趣的是要注意模态形状，尤其是第二个显示出相当大的偏斜（传统高斯混合模型可能在这里失败）。

Mukhopadhyay, S. (2016) 大规模模式识别和数据驱动科学。 https://arxiv.org/abs/1509.06428