为什么 AIC/BIC 标准估计对我的数据的高斯混合密度拟合非常差？

Question

我有一个数据集 (1D) 链接：dataset，其值范围为 21,000 到 8,000,000。当我绘制日志值的直方图时，我可以看到大致有两个峰值。我尝试在 Python 中使用 sklearn 包来拟合高斯混合。我试图根据最低 AIC/BIC 找到最好的 n_components。使用 Full covariance_type，BIC 最好是 44，AIC 是 98（我只测试了 100）。但是一旦我使用这些数字，我的身材就很差。此外，我测试了所有其他 covariance_types，但我无法适应我的数据。我只试了 2 个，我得到了更好的合身度。

这是 44 个组件的图

这是 2 个组件的图

import pandas as pd
import numpy as np
from sklearn.mixture import GaussianMixture as GMM
import matplotlib.pyplot as plt

df = pd.read_excel (r'Data_sets.xlsx',sheet_name="Set1")
b=df['b'].values.reshape(-1,1)
b=np.log(b)
####### finding best n_components ########
k= np.arange(1,100,1)
clfs= [GMM(n,covariance_type='full').fit(b) for n in k]
aics= [clf.aic(b) for clf in clfs]
bics= [clf.bic(b) for clf in clfs]
plt.plot(k,bics,color='orange',marker='.',label='BIC')
plt.plot(k,aics,color='g',label='AIC')
plt.legend()
plt.show()

这是我尝试绘制数据的直方图 + 拟合高斯混合的密度 pdf

clf=GMM(38,covariance_type='full').fit(b)
n, bins, patches = plt.hist(b,bins='auto',density=True,color='#0504aa',alpha=0.7, rwidth=0.85)
xpdf=np.linspace(b.min(),b.max(),len(bins)).reshape(-1,1)
density= np.exp(clf.score_samples(xpdf))
plt.plot(xpdf,density,'-r')
print("Best number of K by BIC is", bics.index(min(bics)))
print("Best number of K by AIC is", aics.index(min(aics)))

这里我绘制了 bins=50 的直方图，顶部直方图用于原始数据集 =3915；根据 BIC 的建议，使用 n_components=44 的 10,000 个样本中的最后一个。看起来 GMM(44) 很合适。

我的问题，导致这些结果的错误在哪里（1）是否因为我的数据不适合高斯混合？ (2) 我的实现错了吗？感谢您提供解决问题的帮助或建议。通过更新（直方图），看起来 GMM 很适合数据。但是，我不明白为什么第一个 plot hist+kde 不合适。我猜是因为 hist 和 kde 都使用不同的 y 比例，但不确定。 谢谢

Answer 1

您的数据似乎由一些非常常见的值和许多相对稀有的值组成。这种区别会混淆混合模型，因此我很想以不同的方式对待它们。如果您不关心这种区别，那么只需像以前一样使用 GMM。我注意到这一点的方法是使用越来越精细的直方图分箱并注意到计数保持不变，表明点质量

我们可以通过numpy.unique 找出这些是什么，例如：

import numpy as np
import pandas as pd

data = pd.read_excel('Data_set.xlsx').values.flatten()

values, counts = np.unique(data, return_counts=True)

# put into a dataframe for nice viewing
uniq = pd.DataFrame(dict(values=values, counts=counts))
uniq.sort_values('counts', ascending=False).head(30)

似乎没有什么好的断点可以使用，所以我任意选择出现超过 10 次的值作为我将特别处理的“流行”值，即作为点质量，我们可以拉出这些做事：

cutoff = 10

popular = set(values[counts > cutoff])
unpopular = [x for x in data if x not in popular]

我们可以将这些绘制在直方图中，并将流行值的叠加计数作为增量峰值给我们：

顶部的箭头表示尖峰在图的顶部（高达 489）在完全主导不受欢迎的值的地方出现，并解释了为什么 GMM 对这些数据如此糟糕地失败（尤其是在对数转换之后）

我将使用高斯 KDE 对“不受欢迎的”数据进行建模，但如果您愿意，也可以使用 GMM。使用 KDE 的一个优点是它是精确的；给定一些数据、内核和带宽，您将始终得到相同的结果。 GMM 要复杂得多，您不太可能每次都得到相同的参数。也就是说，SciPy 中 KDE 的“带宽”参数化是不幸的，但幸运的是我们不需要太多控制，因为分布是最多的

import scipy.stats as sps

kde = sps.gaussian_kde(np.log(unpopular), 0.2)

你可以用这个来说服自己它在做正确的事：

x = np.linspace(11, 16, 501)

plt.hist(np.log(unpopular), 50, density=True)
plt.plot(x, kde(x))

但我不会在此处包含该输出。接下来我们得到一些流行值的汇总统计数据，并定义我们的函数以从中抽取单个样本：

pop_values = values[counts > cutoff]
pop_counts = counts[counts > cutoff]
pop_weights = pop_counts / sum(pop_counts)

pop_prop = sum(pop_counts) / len(data)

def draw_sample():
    if np.random.rand() < pop_prop:
        return np.random.choice(pop_values, p=pop_weights)
    else:
        return int(np.exp(unpop_kde.resample(1)))

samples_10k = [draw_sample() for _ in range(10000)]

最后一行为我们提供了 10k 个样本，我们可以将其绘制成直方图并与原始分布进行比较：

这看起来和我很相似。请注意，这是整体分布，因此以少数点质量为主

如果您只想从“类似的东西”中有效地采样，那么从对数正态基分布中采样的 Dirichlet 过程将起作用并且具有相似的外观属性：

# smaller values will tend to result in more "spikey" classes
alpha = 20

# number of samples to generate
N = 3000

# num components used for finite approximation to dirichlet process, just keep this relatively big
K = 1000

class_values = np.random.lognormal(13, 1, K).astype(int)
class_weight = np.random.dirichlet(np.full(K, alpha/K))
sample_class = np.random.choice(K, N, p=class_weight)
sample_values = class_values[sample_class]

这将从外观相似的分布中生成样本（sample_values 是您可能想要的最终结果），但值会更加多样化。您可能想使用幂律分布来选择样本类（例如 Pitman–Yor 过程）而不是我使用的 Dirichlet，但它们不是内置在 NumPy 中的，因此需要更多代码

Answer 2

发生了一些事情：

1. 根据我的经验，BIC 的额外惩罚通常更好，但使用完整的贝叶斯模型会给您更多控制
2. 您的数据没有很好地被高斯近似。在 40 个组件之后，BIC 不会提高太多，额外的组件只会略微提高拟合度
3. 拟合模型不是确定性的。重新拟合模型，您将导致具有不同 AIC/BIC 的参数略有不同。如果您为每个 n_components 做一些（例如
4. 使用更多 bin 绘制直方图可能有助于在数据中显示更多细节。 bin 的数量取决于数据，但根据我的经验，matplotlib 往往非常保守
5. 当您绘制“直方图”时，我会为xpdf 使用更多数据点（例如 500），这样您就可以更详细地了解正在发生的事情。绘制每个单独的组件可能很有趣/有用。这些可供您使用：

   for i in range(clf.n_components):
     mu_i = clf.means_[i]
     sd_i = np.sqrt(clf.covariances_[i].flatten())
     pr_i = clf.weights_[i]
     density = scipy.stats.norm(mu_i, sd_i).pdf(xpdf) * pr_i
     plt.plot(xpdf, density)
   

   这可以帮助您查看模型将组件放置在何处。对densitys 以上求和与调用clf.score_samples() 相同。下面的演示图：