function [ output_args ] = testGMMandAIC( input_args )
%TESTGMMANDAIC GMMとAICのテスト
% 元のデータの分布がガウス分布より少しでも欠けていたら、AICが推定するカテゴリー数は1にはならないことの確認

%ガウス分布を作る
[Y, Ylabel] = makeDist(1, [-1], [-4], [3000], 'randn', 0);
%ノイズを加える
%noiseNum = 3000;
%noiseDataX = (rand(1, noiseNum) * ( max(Y(1,:)) - min(Y(1,:)) )) + min(Y(1,:));
%noiseDataY = (rand(1, noiseNum) * ( max(Y(2,:)) - min(Y(2,:)) )) + min(Y(2,:));
%noiseData = [noiseDataX ; noiseDataY];
%Y = [Y noiseData];
Y = Y';

%データの右端を少し欠けさせる（SOMの隣接2ノードに対応する入力データはこのような分布をとる）
delnum = find(Y(:,1) > 1.0);
Y(delnum,:) = [];
plot(Y(:,1), Y(:,2), '.');

%データを変換して正規分布に近づける
% まずは負値がなくなるようにずらす
%Yx = Y(:,1) - min(Y(:,1));
%Yy = Y(:,2) - min(Y(:,2));
%Y = [Yx Yy];
%plot(Y(:,1), Y(:,2), '.');
%対数変換
%Y = abs( log(Y) );
%アークサイン変換
%Y = abs( asin(Y) );
%plot(Y(:,1), Y(:,2), '.');

%  GMM
options = statset('Display','final');			% STATS構造体の表示設定
%  最適なカテゴリ数をAICで推定
AIC = zeros(1,10);                               % AICを格納するセル配列
for k = 1:10
    obj_AIC = gmdistribution.fit(Y,k,'Options',options);	% カテゴリ数kで推定
    AIC(k)= obj_AIC.AIC;                        % AICを求める
    % なお、BICの方が多少頑健さが上がる
    %AIC(k)= obj_AIC.BIC;                        % BICを求める
end
AIC
categoryNumAIC = find(AIC == min(AIC));         % AICが最小になるのは何番目かを調べる
disp( strcat( 'カテゴリ数は ', int2str(categoryNumAIC), 'が妥当' ) );

end
