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);

输出如下所示: