# Gaussian Process Mixture

Posted 2014.01.13 01:58

Gaussian Process Mixture Model

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

demo_gpm.m

%% Gaussian process mixture

clc;

clear all;

close all;

%% 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

#### 'Enginius > Machine Learning' 카테고리의 다른 글

 Reinforcement Learning  (5) 2014.02.07 2014.01.23 2014.01.13 2014.01.09 2013.10.15 2013.07.22
« PREV : 1 : ··· : 75 : 76 : 77 : 78 : 79 : 80 : 81 : 82 : 83 : ··· : 142 : NEXT »