Gaussian Mixture Regression-Code
Zhao Zhao
  2023-05-13 17:21:06 2023-05-13 17:21    

Handbook of Robotics Chapter 59: Robot Programming by Demonstration

A dataset is defined by observations of sensory data changing through time (e.g., joint angle trajectories, hand paths), where each datapoint consists of a temporal value and a spatial vector . The dataset is modelled by a Gaussian Mixture Model (GMM) of components, defined by the probability density function where are prior probabilities and are Gaussian distributions defined by mean vectors and covariance matrices , whose temporal and spatial components can be represented separately as For each component , the expected distribution of given the temporal value is defined by By considering the complete GMM, the expected distribution is defined by where is the probability of the component to be responsible for , i.e., By using the linear transformation properties of Gaussian distributions, an estimation of the conditional expectation of give is thus defined by , where the parameters of the Gaussian distribution are defined by By evaluating at different time steps , a generalized form of the motions and associated covariance matrices describing the constraints are computed.

Mixture models for the analysis, edition, and synthesis of continuous time series

The proof is in A tutorial on task-parameterized movement learning and retrieval

In this work the GMM is given by

Method of Moments

https://www.math.arizona.edu/~jwatkins/M_moments.pdf

https://web.stanford.edu/class/archive/cs/cs109/cs109.1218/files/student_drive/7.3.pdf

Gaussian mixture regression and classification https://scholarship.rice.edu/bitstream/handle/1911/18710/3122551.PDF?sequence=1

Here, we derive the method of moments estimations of a 2-component GMM. Given the pdf of Compute the mean and variance of .

First, derive and : 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
clear
clc
close all

nbStates = 4;  % 高斯分布个数
load('data1.mat'); %load 'Data'
% data : 3次示教数据，每次100个点，data(1, : )--time index, data(2, : )--x1 data, data(3, : )--x2 data
nbVar = size(Data,1);   % 3个变量，第一个是输入，后两个是输出

[Priors, Mu, Sigma] = EM_init_kmeans(Data, nbStates);
[Priors, Mu, Sigma] = EM(Data, Priors, Mu, Sigma);

% Outputs ----------------------------------------------------------------
%   o Priors:   1 x K array representing the prior probabilities of the
%               K GMM components.
%   o Mu:       D x K array representing the centers of the K GMM components.
%   o Sigma:    D x D x K array representing the covariance matrices of the 
%               K GMM components.


expData(1,:) = linspace(min(Data(1,:)), max(Data(1,:)), 100);   % 输入

x = expData(1,:) ; % input
in = 1; % 1 x P array representing the dimensions to consider as  inputs
out = [2:nbVar];  % 输出的数据维度   1 x Q array representing the dimensions to consider as
%              outputs (D=P+Q).


nbData = length(x);  % 输入点数
nbVar = size(Mu,1);  %  3个变量，第一个是输入，后两个是输出
nbStates = size(Sigma,3);  % K = 4


%%  Compute the influence of each GMM component, given input x

for k = 1:nbStates
        Pxi(k, :) = Priors(k).* mvnpdf(x', Mu(in, k), Sigma(in ,in, k));
%         Pxi(:,k) = Priors(k).*gaussPDF(x, Mu(in,k), Sigma(in,in,k));
end
r = Pxi./sum(Pxi);  % 这里sum(r)是列相加 % K x N
% r = (Pxi./repmat(sum(Pxi,2)+realmin,1,nbStates))';
%% Compute expected output distribution, given input x
Sigma_y = zeros(length(out), length(out), nbData);    % 2 x 2x 4
y = zeros(length(out), nbData);
for n = 1:nbData
    % Compute expected means y, given input x
    for k = 1:nbStates
         yj_tmp =  Mu(out,k) + Sigma(out,in,k)*inv(Sigma(in,in,k)) * (x(:,n)-Mu(in,k));
         y(:,n) = y(:,n) + r(k, n).*yj_tmp; 
    end
  % Compute expected covariance matrices Sigma_y, given input x
%     for k = 1:nbStates
%          Sigmaj_y_tmp = Sigma(out,out,k) - (Sigma(out,in,k)*inv(Sigma(in,in,k))*Sigma(in,out,k));
%          Sigma_y(:,:,n) = Sigma_y(:,:,n) + r(k, n)^2.* Sigmaj_y_tmp;
%     end
    for k = 1:nbStates
        yj_tmp =  Mu(out,k) + Sigma(out,in,k)*inv(Sigma(in,in,k)) * (x(:,n)-Mu(in,k));
        Sigmaj_y_tmp = Sigma(out,out,k) - (Sigma(out,in,k)*inv(Sigma(in,in,k))*Sigma(in,out,k));
        Sigma_y(:,:,n) = Sigma_y(:,:,n) + r(k, n).* (Sigmaj_y_tmp+yj_tmp*yj_tmp');
    end
     Sigma_y(:,:,n) = Sigma_y(:,:,n)- y(:,n)*y(:,n)';
end

expData(2:nbVar,:) = y;
expSigma = Sigma_y ;
% [expData(2:nbVar,:), expSigma] = GMR(Priors, Mu, Sigma,  expData(1,:), [1], [2:nbVar]);

%% GMR%% Plot of the data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure('position',[10,10,1000,800],'name','GMM-GMR-demo1');
%plot 1D
for n=1:nbVar-1
  subplot(3*(nbVar-1),2,(n-1)*2+1); hold on;
  plot(Data(1,:), Data(n+1,:), 'x', 'markerSize', 4, 'color', [.3 .3 .3]);
  axis([min(Data(1,:)) max(Data(1,:)) min(Data(n+1,:))-0.01 max(Data(n+1,:))+0.01]);
  xlabel('t','fontsize',16); ylabel(['x_' num2str(n)],'fontsize',16);
end
%plot 2D
subplot(3*(nbVar-1),2,[2:2:2*(nbVar-1)]); hold on;
plot(Data(2,:), Data(3,:), 'x', 'markerSize', 4, 'color', [.3 .3 .3]);
axis([min(Data(2,:))-0.01 max(Data(2,:))+0.01 min(Data(3,:))-0.01 max(Data(3,:))+0.01]);
xlabel('x_1','fontsize',16); ylabel('x_2','fontsize',16);

%% Plot of the GMM encoding results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%plot 1D
for n=1:nbVar-1
  subplot(3*(nbVar-1),2,4+(n-1)*2+1); hold on;
  plotGMM(Mu([1,n+1],:), Sigma([1,n+1],[1,n+1],:), [0 .8 0], 1);
  axis([min(Data(1,:)) max(Data(1,:)) min(Data(n+1,:))-0.01 max(Data(n+1,:))+0.01]);
  xlabel('t','fontsize',16); ylabel(['x_' num2str(n)],'fontsize',16);
end
%plot 2D
subplot(3*(nbVar-1),2,4+[2:2:2*(nbVar-1)]); hold on;
plotGMM(Mu([2,3],:), Sigma([2,3],[2,3],:), [0 .8 0], 1);
axis([min(Data(2,:))-0.01 max(Data(2,:))+0.01 min(Data(3,:))-0.01 max(Data(3,:))+0.01]);
xlabel('x_1','fontsize',16); ylabel('x_2','fontsize',16);

%% Plot of the GMR regression results
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%plot 1D
for n=1:nbVar-1
  subplot(3*(nbVar-1),2,8+(n-1)*2+1); hold on;
  plotGMM(expData([1,n+1],:), expSigma(n,n,:), [0 0 .8], 3);
  axis([min(Data(1,:)) max(Data(1,:)) min(Data(n+1,:))-0.01 max(Data(n+1,:))+0.01]);
  xlabel('t','fontsize',16); ylabel(['x_' num2str(n)],'fontsize',16);
end
%plot 2D
subplot(3*(nbVar-1),2,8+[2:2:2*(nbVar-1)]); hold on;
plotGMM(expData([2,3],:), expSigma([1,2],[1,2],:), [0 0 .8], 2);
axis([min(Data(2,:))-0.01 max(Data(2,:))+0.01 min(Data(3,:))-0.01 max(Data(3,:))+0.01]);
xlabel('x_1','fontsize',16); ylabel('x_2','fontsize',16);

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
function plotGMM(Mu, Sigma, color, display_mode);
%
% This function plots a representation of the components (means and 
% covariance amtrices) of a Gaussian Mixture Model (GMM) or a
% Gaussian Mixture Regression (GMR).
%
% Author:	Sylvain Calinon, 2009
%			http://programming-by-demonstration.org
%
% Inputs -----------------------------------------------------------------
%   o Mu:           D x K array representing the centers of the K GMM components.
%   o Sigma:        D x D x K array representing the covariance matrices of the 
%                   K GMM components.
%   o color:        Color used for the representation
%   o display_mode: Display mode (1 is used for a GMM, 2 is used for a GMR
%                   with a 2D representation and 3 is used for a GMR with a
%                   1D representation).

nbData = size(Mu,2);
lightcolor = color + [0.6,0.6,0.6];
lightcolor(find(lightcolor>1.0)) = 1.0;

if display_mode==1
  nbDrawingSeg = 40;
  t = linspace(-pi, pi, nbDrawingSeg)';
  for j=1:nbData
    stdev = sqrtm(3.0.*Sigma(:,:,j));
    X = [cos(t) sin(t)] * real(stdev) + repmat(Mu(:,j)',nbDrawingSeg,1);
    patch(X(:,1), X(:,2), lightcolor, 'lineWidth', 2, 'EdgeColor', color);
    plot(Mu(1,:), Mu(2,:), 'x', 'lineWidth', 2, 'color', color);
  end
elseif display_mode==2
  nbDrawingSeg = 40;
  t = linspace(-pi, pi, nbDrawingSeg)';
  for j=1:nbData
    stdev = sqrtm(3.0.*Sigma(:,:,j));
    X = [cos(t) sin(t)] * real(stdev) + repmat(Mu(:,j)',nbDrawingSeg,1);
    patch(X(:,1), X(:,2), lightcolor, 'LineStyle', 'none');
  end
  plot(Mu(1,:), Mu(2,:), '-', 'lineWidth', 3, 'color', color);
elseif display_mode==3
  for j=1:nbData
    ymax(j) = Mu(2,j) + sqrtm(3.*Sigma(1,1,j));
    ymin(j) = Mu(2,j) - sqrtm(3.*Sigma(1,1,j));
  end
  patch([Mu(1,1:end) Mu(1,end:-1:1)], [ymax(1:end) ymin(end:-1:1)], lightcolor, 'LineStyle', 'none');
  plot(Mu(1,:), Mu(2,:), '-', 'lineWidth', 3, 'color', color); 
end