Gaussian Process Mixture

Gaussian Process Mixture Model 

물론 아래 프로그램을 돌리기 위해선 GPRobotLib이 필요하다. 

demo_gpm.m

%% Gaussian process mixture

clc;

clear all;

close all;

addpath('../CustomLibraries/GPRobotLib');

%% PLOT OPTION

plot_ref   = 1;

plot_gpr   = 1;

plot_gpm   = 1;

%% Generate test data 

refX = linspace(0, 50, 200);

refY = f(refX);

gpRes = 50;

gpX = linspace(0, 50, gpRes);

gpY = f(gpX);

if plot_ref

    figure(1);

    hold on;

    href = plot(refX, refY, 'b-' ...

        , 'LineWidth', 2);

    hgp = plot(gpX, gpY, 'ro' ...

        , 'LineWidth', 2);

    hold off;

    axis([0 50 -1.5 1.5]);

    legend([href, hgp], 'Ref', 'GP data');

    title('Generated data', 'FontSize', 15);

    xlabel('X', 'FontSize', 15);

    ylabel('Y', 'FontSize', 15);

end

%% GPR을 한다. 

gp_data.nr_data = gpRes;

gp_data.input  = gpX;

gp_data.output = gpY;

gp_data.Z = linspace(0, 50, 1000)';

hyp.sig2f  = 1E+0;

hyp.sig2w  = 1E-4;

hyp.sig2x  = 5E+0;

[gp_struct]...

    = initGP1D(gp_data ...

    , @gp_kernel_1D_se, hyp);

gp_result = gpr(gp_struct, gp_data.Z);

if plot_gpr

    hfig = figure(2);

    set(hfig, 'Position', [100 100 900 700]);

    hold on;

    href = plot(gp_data.input, gp_data.output, 'ro-' ...

        , 'LineWidth', 2);

    hgp = plot(gp_data.Z, gp_result.mean, 'b-' ...

        , 'LineWidth', 2);

    v_gain = 1E3;

    plot(gp_data.Z, gp_result.mean + v_gain*gp_result.var, 'g--' ...

        , 'LineWidth', 1);

    plot(gp_data.Z, gp_result.mean - v_gain*gp_result.var, 'g--' ...

        , 'LineWidth', 1);

    hold off; 

    legend([href, hgp], 'Reference', 'GPR');

    axis([0 50 -1.5 1.5]);

    title('Gaussian Process Regression Result', 'FontSize', 15);

    xlabel('X', 'FontSize', 15);

    ylabel('Y', 'FontSize', 15);

end

%% Do Gaussian Process Mixture (GPM)

% Init Gaussian Process Mixture (GPM)

gpm_struct = initGPM(gp_data, @gp_kernel_1D_se, hyp);

if plot_gpm

    fimgpe = figure(3);

    set(fimgpe, 'Position', [100 100 900 700]);

    img_cnt = 0;

    colors = varyColor(gpm_struct.nr_cluster);

end

% Loop

loop_struct = initLoop(1E2);

while loop_struct.flag 

    % Terminate condition 

    loop_struct.tick = loop_struct.tick + 1;

    if loop_struct.tick > loop_struct.maxTick

        break; 

    end 

    fprintf('[%d/%d] \n', loop_struct.tick, loop_struct.maxTick);

    

    % 그림을 그려서 확인한다. 

    if plot_gpm

        clf;

        hold on;

        for i = 1:gpm_struct.nr_data

            curr_probs = gpm_struct.class_prob(i, :);

            curr_color = zeros(1, 3);

            for j = 1:gpm_struct.nr_cluster

                curr_color = curr_color + curr_probs(j)*colors(j, :);

            end

            % color range를 0~1로 한다. 

            curr_color(find(curr_color<0)) = 0;

            curr_color(find(curr_color>1)) = 1;

            plot(gpm_struct.input(i), gpm_struct.output(i), 'o' ...

                , 'MarkerSize', 10 ...

                , 'MarkerFaceColor', curr_color);

        end

        hold off;

        title( sprintf('[%d/%d] ' ...

            , loop_struct.tick, loop_struct.maxTick) ...

            , 'FontSize', 15);

        drawnow;

        % pause();

    end

    

    % 모든 데이터에 대해서 각 cluster의 확률을 update한다. 

    class_prob_temp = gpm_struct.class_prob;

    % 각 cluster에 해당하는 데이터를 모은다. 

    [~, max_indices] = max(class_prob_temp');

    

    for i = 1:gpm_struct.nr_data

        % 각 데이터에 대해서 

        x_i = gpm_struct.input(i, :);

        y_i = gpm_struct.output(i, :);

        % 현재 데이터에 대한 각 class의 확률 

        class_prob_temp_i = class_prob_temp(i, :);

        

        % 각 cluster를 구한다. 

        gp_clusters_i = cell(gpm_struct.nr_cluster, 1);

        for j = 1:gpm_struct.nr_cluster 

            cluster_indices = find(max_indices == j);

            % cluster_indices 에서 i를 제외한다. 

            cluster_indices(find(cluster_indices == i)) = [];

            

            % GP 데이터를 만든다. 

            gp_data_temp.input = gpm_struct.input(cluster_indices);

            gp_data_temp.output = gpm_struct.output(cluster_indices);

            gp_data_temp.nr_data = length(cluster_indices);

            hyp_temp = hyp;        

            gp_clusters_i{j} = initGP1D(gp_data_temp, @gp_kernel_1D_se, hyp);

        end

        

        % i번째 data의 class prob을 구한다. 

        class_prob_i = zeros(1, gpm_struct.nr_cluster);

        for j = 1:gpm_struct.nr_cluster

            if gp_clusters_i{j}.nr_data == 0

                % cluster 

                class_prob_i(j) = 0;

            else

                gp_result = gpr(gp_clusters_i{j}, x_i);

                mean_i_j = gp_result.mean;

                var_i_j = gp_result.var;

                class_prob_i(j) = class_prob_temp_i(j)*mvnpdf(y_i, mean_i_j, var_i_j);

            end

        end

        class_prob_i = class_prob_i / sum(class_prob_i);

        

        % Update한다. 

        gpm_struct.class_prob(i, :) = class_prob_i;

    end

end

%%

%%

fprintf('Thats all folks \n');

initGPM.m

function gpm_struct = initGPM(gp_data, gp_kernel, hyp)

% initialize imgpe struct

gpm_struct.nr_data = gp_data.nr_data; 

gpm_struct.input   = reshape(gp_data.input, gpm_struct.nr_data, []); 

gpm_struct.output  = reshape(gp_data.output, gpm_struct.nr_data, []); 

% 현재 GPM의 cluster의 수 (k)

gpm_struct.nr_cluster = 5;

% 각 데이터가 소속된 cluster (c_i)

% 일단 1번 cluster에 모든 데이터가 들어가있게 한다. 

gpm_struct.class_prob = ones(gp_data.nr_data, gpm_struct.nr_cluster);

for i = 1:gp_data.nr_data

    curr_rand = rand(1, gpm_struct.nr_cluster);

    curr_rand = curr_rand / sum(curr_rand);

    gpm_struct.class_prob(i, :) = curr_rand;

end

% Concentration parameter

gpm_struct.alpha = 1E-1;

gpm_struct.gkernel = @(X1, X2)gp_kernel(X1, X2, hyp);

gpm_struct.rawkernel = gp_kernel;

f.m

function ret = f(X)

nrX = length(X);

ret = zeros(nrX, 1);

for i = 1:nrX

    if X(i) < 10

        ret(i) = 1;

    elseif X(i) < 20

        ret(i) = -1;

    elseif X(i) < 30

        ret(i) = 1;

    elseif X(i) < 40

        ret(i) = -1;

    else

        ret(i) = 1;

    end

end

multinomialRand.m

function [ret_mat, ret_vec]...

    = multinomialRand(n, p)

%

% [a, b] = multinomialRand(1, [0.1, 0.2, 0.7])

% 결과가 다음과 같이 나온다. 

% a = [0 0 1], b = 3

% a = [0 0 1], b = 3

% a = [1 0 0], b = 1

% a = [0 1 0], b = 2

% a = [0 0 1], b = 3

%

% 혹시 모르니 p를 normalize한다. 

p = p/sum(p);

nrVar = length(p);

ret_mat = zeros(n, nrVar);

ret_vec = zeros(n);

for i = 1:n

    idx = find(mnrnd(1, p), 1);

    ret_mat(i, idx) = 1;

    ret_vec(i) = idx;

end