Matlab：使用许多高斯和 AIC 评估拟合直方图数据

Question

考虑这个代码示例，根据 Akaike 准则从改变拟合高斯数的数据中获得最佳拟合

MU1 = [1];
SIGMA1 = [2];
MU2 = [-3];
SIGMA2 = [1 ];
X = [mvnrnd(MU1,SIGMA1,1000);mvnrnd(MU2,SIGMA2,1000)];
AIC = zeros(1,4);
obj = cell(1,4);
options = statset('Display','final');
for k = 1:4
 obj{k} = gmdistribution.fit(X,k,'Options',options);
 AIC(k)= obj{k}.AIC;
end
[minAIC,numComponents] = min(AIC)

我想做同样的事情，但数据是以直方图的形式给出的（例如数据http://pastebin.com/embed_js.php?i=1mNRuEHZ）。

在这种情况下，在matlab中实现相同程序的最直接方法是什么？

Answer 1

如果我说得对，那么您的问题是在已编译为直方图的数据（因此观察的数量与观察的实际值配对）和原始的个人观察之间进行转换。当然，在编译直方图的时候，你丢失了两件事：

1. 订购。您不知道原始数据中观察的顺序是什么，这可能并不重要，前提是您的观察是独立的。另外，我得到 gmdistribution.fit() 的方式无论如何都没有考虑顺序。

2. 分辨率。当您创建直方图时，您需要对数据进行分箱，这会使您失去精度，可以这么说，因为无法从分箱中恢复您的观察值的精确值。

一旦您意识到这一点，您就可以从您的直方图数据中创建一个“观察向量”。说，X1 是您的直方图数据（Nx2 向量）。如果你这样做了

invX = cell2mat(arrayfun(@(x,y) repmat(y,1,x), abs(int16(1000*X1(:, 2)))', X1(:, 1)', ...
    'UniformOutput', false))';

您会得到一个包含单个观察值的向量，就像您示例中的 X 一样。

请注意，您必须先将 bin 计数转换为整数。在这一步，由于给定数据的精度很高，我不得不四舍五入以使我的机器可以计算。然而，最终的结果似乎相当合理。

另外请注意，我使用了绝对值，在您的直方图数据中有些情况是您的数据实际上是负数，这对于直方图显然没有意义。

最后但并非最不重要的一点是，您必须将拟合过程的迭代次数更改为 1000。生成下图的最终代码为

MU1 = [1];
SIGMA1 = [2];
MU2 = [-3];
SIGMA2 = [1 ];
X = [mvnrnd(MU1,SIGMA1,1000);mvnrnd(MU2,SIGMA2,1000)];
X = X1(:, 2);
invX = cell2mat(arrayfun(@(x,y) repmat(y,1,x), abs(int16(1000*X1(:, 2)))', X1(:, 1)', ...
    'UniformOutput', false))'; %'
X = invX;
AIC = zeros(1,4);
obj = cell(1,4);
options = statset('Display','final', 'MaxIter', 1000);
for k = 1:4
 obj{k} = gmdistribution.fit(X,k,'Options',options);
 AIC(k)= obj{k}.AIC;
end
[minAIC,numComponents] = min(AIC);

hold on;
plot(linspace(-1, 2, length(X1(:, 2))), abs(X1(:, 2)), 'LineWidth', 2)
plot(x, pd/max(pd)*double(max(abs(X1(:, 2)))), 'LineWidth', 5);
h = legend('Original data', 'PDF');
set(h,'FontSize',32);

输出如下所示：