Skip to content
Projects
Groups
Snippets
Help
Loading...
Help
Contribute to GitLab
Sign in
Toggle navigation
pbdlibmatlab
Project
Project
Details
Activity
Cycle Analytics
Repository
Repository
Files
Commits
Branches
Tags
Contributors
Graph
Compare
Charts
Issues
0
Issues
0
List
Board
Labels
Milestones
Merge Requests
0
Merge Requests
0
CI / CD
CI / CD
Pipelines
Jobs
Schedules
Charts
Members
Members
Collapse sidebar
Close sidebar
Activity
Graph
Charts
Create a new issue
Jobs
Commits
Issue Boards
Open sidebar
rli
pbdlibmatlab
Commits
2cb8784d
Commit
2cb8784d
authored
Sep 16, 2015
by
Sylvain Calinon
Browse files
Options
Browse Files
Download
Email Patches
Plain Diff
Homogenization of headers in m_fcts folder
parent
b738860a
Hide whitespace changes
Inline
Sidebyside
Showing
38 changed files
with
1316 additions
and
414 deletions
+1316
414
demo_stdPGMM01.m
demo_stdPGMM01.m
+13
2
DTW.m
m_fcts/DTW.m
+42
18
EM_GMM.m
m_fcts/EM_GMM.m
+35
10
EM_HDGMM.m
m_fcts/EM_HDGMM.m
+48
9
EM_MFA.m
m_fcts/EM_MFA.m
+42
13
EM_MPPCA.m
m_fcts/EM_MPPCA.m
+41
12
EM_TPGMM.m
m_fcts/EM_TPGMM.m
+33
23
EM_stdPGMM.m
m_fcts/EM_stdPGMM.m
+33
25
EM_tensorGMM.m
m_fcts/EM_tensorGMM.m
+33
23
EM_tensorHDGMM.m
m_fcts/EM_tensorHDGMM.m
+39
14
EM_tensorMFA.m
m_fcts/EM_tensorMFA.m
+42
16
EM_tensorMPPCA.m
m_fcts/EM_tensorMPPCA.m
+43
18
GMR.m
m_fcts/GMR.m
+31
21
GPR.m
m_fcts/GPR.m
+32
5
constructPHI.m
m_fcts/constructPHI.m
+36
9
estimateAttractorPath.m
m_fcts/estimateAttractorPath.m
+29
16
gaussPDF.m
m_fcts/gaussPDF.m
+29
4
init_GMM_kmeans.m
m_fcts/init_GMM_kmeans.m
+31
12
init_GMM_timeBased.m
m_fcts/init_GMM_timeBased.m
+33
15
init_TPGMM_kmeans.m
m_fcts/init_TPGMM_kmeans.m
+35
12
init_TPGMM_timeBased.m
m_fcts/init_TPGMM_timeBased.m
+32
4
init_tensorGMM_kmeans.m
m_fcts/init_tensorGMM_kmeans.m
+35
12
init_tensorGMM_timeBased.m
m_fcts/init_tensorGMM_timeBased.m
+32
4
kmeansClustering.m
m_fcts/kmeansClustering.m
+31
6
plot2DArrow.m
m_fcts/plot2DArrow.m
+32
2
plotArm.m
m_fcts/plotArm.m
+31
1
plotArmBasis.m
m_fcts/plotArmBasis.m
+33
3
plotArmLink.m
m_fcts/plotArmLink.m
+34
5
plotBimanualRobot.m
m_fcts/plotBimanualRobot.m
+45
16
plotGMM.m
m_fcts/plotGMM.m
+29
2
productTPGMM.m
m_fcts/productTPGMM.m
+34
6
productTPGMM0.m
m_fcts/productTPGMM0.m
+32
3
reproduction_DS.m
m_fcts/reproduction_DS.m
+31
10
reproduction_LQR_finiteHorizon.m
m_fcts/reproduction_LQR_finiteHorizon.m
+38
15
reproduction_LQR_infiniteHorizon.m
m_fcts/reproduction_LQR_infiniteHorizon.m
+35
13
reproduction_TPGMM.m
m_fcts/reproduction_TPGMM.m
+28
16
solveAlgebraicRiccati_Schur.m
m_fcts/solveAlgebraicRiccati_Schur.m
+46
13
solveAlgebraicRiccati_eig.m
m_fcts/solveAlgebraicRiccati_eig.m
+38
6
No files found.
demo_stdPGMM01.m
View file @
2cb8784d
function
demo_stdPGMM01
% Parametric Gaussian mixture model (PGMM) used for task adaptation,
% with DSGMR employed to retrieve continuous movements.
% Parametric Gaussian mixture model (PGMM) used for task adaptation, with DSGMR employed
% to retrieve continuous movements. The approach is inspired by Wilson and Bobick (1999),
% with an implementation applied to the special case of Gaussian mixture models (GMM).
%
% Writing code takes time. Polishing it and making it available to others takes longer!
% If some parts of the code were useful for your research of for a better understanding
...
...
@@ 13,6 +14,16 @@ function demo_stdPGMM01
% journal="Intelligent Service Robotics",
% year="2015"
% }
%
% @article{Wilson99,
% author="Wilson, A. D. and Bobick, A. F.",
% title="Parametric Hidden {M}arkov Models for Gesture Recognition",
% journal="{IEEE} Trans. on Pattern Analysis and Machine Intelligence",
% year="1999",
% volume="21",
% number="9",
% pages="884900"
% }
%
% Copyright (c) 2015 Idiap Research Institute, http://idiap.ch/
% Written by Sylvain Calinon, http://calinon.ch/
...
...
m_fcts/DTW.m
View file @
2cb8784d
function
[
x_new
,
y_new
,
p
]
=
DTW
(
x
,
y
,
w
)
%Trajectory realignment through dynamic time warping
%Sylvain Calinon, 2015
function
[
x2
,
y2
,
p
]
=
DTW
(
x
,
y
,
w
)
% Trajectory realignment through dynamic time warping.
%
% Writing code takes time. Polishing it and making it available to others takes longer!
% If some parts of the code were useful for your research of for a better understanding
% of the algorithms, please reward the authors by citing the related publications,
% and consider making your own research available in this way.
%
% @article{Calinon15,
% author="Calinon, S.",
% title="A Tutorial on TaskParameterized Movement Learning and Retrieval",
% journal="Intelligent Service Robotics",
% year="2015"
% }
%
% Copyright (c) 2015 Idiap Research Institute, http://idiap.ch/
% Written by Sylvain Calinon, http://calinon.ch/
%
% This file is part of PbDlib, http://www.idiap.ch/software/pbdlib/
%
% PbDlib is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License version 3 as
% published by the Free Software Foundation.
%
% PbDlib is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with PbDlib. If not, see <http://www.gnu.org/licenses/>.
if
nargin
<
3
w
=
Inf
;
end
n
x
=
size
(
x
,
2
);
n
y
=
size
(
y
,
2
);
s
x
=
size
(
x
,
2
);
s
y
=
size
(
y
,
2
);
w
=
max
(
w
,
abs
(
nx

n
y
));
w
=
max
(
w
,
abs
(
sx

s
y
));
%Initialization
D
=
ones
(
nx
+
1
,
n
y
+
1
)
*
Inf
;
D
=
ones
(
sx
+
1
,
s
y
+
1
)
*
Inf
;
D
(
1
,
1
)
=
0
;
%DP loop
for
i
=
1
:
nx
for
j
=
max
(
i

w
,
1
):
min
(
i
+
w
,
n
y
)
for
j
=
max
(
i

w
,
1
):
min
(
i
+
w
,
s
y
)
D
(
i
+
1
,
j
+
1
)
=
norm
(
x
(:,
i
)

y
(:,
j
))
+
min
([
D
(
i
,
j
+
1
),
D
(
i
+
1
,
j
),
D
(
i
,
j
)]);
end
end
i
=
nx
+
1
;
j
=
n
y
+
1
;
i
=
sx
+
1
;
j
=
s
y
+
1
;
p
=
[];
while
i
>
1
&&
j
>
1
[
~
,
id
]
=
min
([
D
(
i
,
j

1
),
D
(
i

1
,
j
),
D
(
i

1
,
j

1
)]);
...
...
@@ 39,14 +68,9 @@ end
p
=
fliplr
(
p
(:,
1
:
end

1
)

1
);
x
_new
=
x
(:,
p
(
1
,:));
y
_new
=
y
(:,
p
(
2
,:));
x
2
=
x
(:,
p
(
1
,:));
y
2
=
y
(:,
p
(
2
,:));
%Resampling
x_new
=
spline
(
1
:
size
(
x_new
,
2
),
x_new
,
linspace
(
1
,
size
(
x_new
,
2
),
nx
));
y_new
=
spline
(
1
:
size
(
y_new
,
2
),
y_new
,
linspace
(
1
,
size
(
y_new
,
2
),
nx
));
x2
=
spline
(
1
:
size
(
x2
,
2
),
x2
,
linspace
(
1
,
size
(
x2
,
2
),
sx
));
y2
=
spline
(
1
:
size
(
y2
,
2
),
y2
,
linspace
(
1
,
size
(
y2
,
2
),
sx
));
m_fcts/EM_GMM.m
View file @
2cb8784d
function
[
model
,
GAMMA2
,
LL
]
=
EM_GMM
(
Data
,
model
)
% Training of a Gaussian mixture model (GMM) with an expectationmaximization (EM) algorithm.
%
% Author: Sylvain Calinon, 2014
% http://programmingbydemonstration.org/SylvainCalinon
% Writing code takes time. Polishing it and making it available to others takes longer!
% If some parts of the code were useful for your research of for a better understanding
% of the algorithms, please reward the authors by citing the related publications,
% and consider making your own research available in this way.
%
% @article{Calinon15,
% author="Calinon, S.",
% title="A Tutorial on TaskParameterized Movement Learning and Retrieval",
% journal="Intelligent Service Robotics",
% year="2015"
% }
%
% Copyright (c) 2015 Idiap Research Institute, http://idiap.ch/
% Written by Sylvain Calinon, http://calinon.ch/
%
% This file is part of PbDlib, http://www.idiap.ch/software/pbdlib/
%
% PbDlib is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License version 3 as
% published by the Free Software Foundation.
%
% PbDlib is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with PbDlib. If not, see <http://www.gnu.org/licenses/>.
%Parameters of the EM algorithm
nbMinSteps
=
5
;
%Minimum number of iterations allowed
...
...
@@ 10,8 +37,8 @@ nbMaxSteps = 100; %Maximum number of iterations allowed
maxDiffLL
=
1E4
;
%Likelihood increase threshold to stop the algorithm
nbData
=
size
(
Data
,
2
);
%diagRegularizationFactor = 1E6; %Regularization term is optional
, see Eq. (2.1.2) in doc/TechnicalReport.pdf
diagRegularizationFactor
=
1E4
;
%Regularization term is optional
, see Eq. (2.1.2) in doc/TechnicalReport.pdf
%diagRegularizationFactor = 1E6; %Regularization term is optional
diagRegularizationFactor
=
1E4
;
%Regularization term is optional
for
nbIter
=
1
:
nbMaxSteps
fprintf
(
'.'
);
...
...
@@ 22,13 +49,13 @@ for nbIter=1:nbMaxSteps
%Mstep
for
i
=
1
:
model
.
nbStates
%Update Priors
, see Eq. (2.0.6) in doc/TechnicalReport.pdf
%Update Priors
model
.
Priors
(
i
)
=
sum
(
GAMMA
(
i
,:))
/
nbData
;
%Update Mu
, see Eq. (2.0.7) in doc/TechnicalReport.pdf
%Update Mu
model
.
Mu
(:,
i
)
=
Data
*
GAMMA2
(
i
,:)
'
;
%Update Sigma
, see Eq. (2.0.8) in doc/TechnicalReport.pdf (regularization term is optional, see Eq. (2.1.2))
%Update Sigma
DataTmp
=
Data

repmat
(
model
.
Mu
(:,
i
),
1
,
nbData
);
model
.
Sigma
(:,:,
i
)
=
DataTmp
*
diag
(
GAMMA2
(
i
,:))
*
DataTmp
'
+
eye
(
size
(
Data
,
1
))
*
diagRegularizationFactor
;
end
...
...
@@ 46,9 +73,9 @@ end
disp
([
'The maximum number of '
num2str
(
nbMaxSteps
)
' EM iterations has been reached.'
]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function
[
L
,
GAMMA
]
=
computeGamma
(
Data
,
model
)
%See Eq. (2.0.5) in doc/TechnicalReport.pdf
L
=
zeros
(
model
.
nbStates
,
size
(
Data
,
2
));
for
i
=
1
:
model
.
nbStates
L
(
i
,:)
=
model
.
Priors
(
i
)
*
gaussPDF
(
Data
,
model
.
Mu
(:,
i
),
model
.
Sigma
(:,:,
i
));
...
...
@@ 56,5 +83,3 @@ end
GAMMA
=
L
.
/
repmat
(
sum
(
L
,
1
)
+
realmin
,
model
.
nbStates
,
1
);
end
m_fcts/EM_HDGMM.m
View file @
2cb8784d
function
[
model
,
GAMMA2
]
=
EM_HDGMM
(
Data
,
model
)
%EM for High Dimensional Data Clustering (HDDC, HDGMM) model proposed by Bouveyron (2007)
%Sylvain Calinon, 2015
% EM for High Dimensional Data Clustering (HDDC, HDGMM) model proposed by Bouveyron (2007).
%
% Writing code takes time. Polishing it and making it available to others takes longer!
% If some parts of the code were useful for your research of for a better understanding
% of the algorithms, please reward the authors by citing the related publications,
% and consider making your own research available in this way.
%
% @article{Calinon15,
% author="Calinon, S.",
% title="A Tutorial on TaskParameterized Movement Learning and Retrieval",
% journal="Intelligent Service Robotics",
% year="2015"
% }
%
% @article{Bouveyron07,
% author = "Bouveyron, C. and Girard, S. and Schmid, C.",
% title = "Highdimensional data clustering",
% journal = "Computational Statistics and Data Analysis",
% year = "2007",
% volume = "52",
% number = "1",
% pages = "502519"
% }
%
% Copyright (c) 2015 Idiap Research Institute, http://idiap.ch/
% Written by Sylvain Calinon, http://calinon.ch/
%
% This file is part of PbDlib, http://www.idiap.ch/software/pbdlib/
%
% PbDlib is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License version 3 as
% published by the Free Software Foundation.
%
% PbDlib is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with PbDlib. If not, see <http://www.gnu.org/licenses/>.
%Parameters of the EM iterations
nbMinSteps
=
5
;
%Minimum number of iterations allowed
...
...
@@ 8,7 +47,7 @@ nbMaxSteps = 100; %Maximum number of iterations allowed
maxDiffLL
=
1E4
;
%Likelihood increase threshold to stop the algorithm
nbData
=
size
(
Data
,
2
);
diagRegularizationFactor
=
1E8
;
%Regularization term is optional
, see Eq. (2.1.2) in doc/TechnicalReport.pdf
diagRegularizationFactor
=
1E8
;
%Regularization term is optional
%EM loop
for
nbIter
=
1
:
nbMaxSteps
...
...
@@ 19,19 +58,19 @@ for nbIter=1:nbMaxSteps
GAMMA2
=
GAMMA
.
/
repmat
(
sum
(
GAMMA
,
2
),
1
,
nbData
);
%Mstep
%Update Priors
, see Eq. (2.0.6) in doc/TechnicalReport.pdf
%Update Priors
model
.
Priors
=
sum
(
GAMMA
,
2
)/
nbData
;
%Update Mu
, see Eq. (2.0.7) in doc/TechnicalReport.pdf
%Update Mu
model
.
Mu
=
Data
*
GAMMA2
'
;
%Update factor analyser params
%Update factor analyser param
eter
s
for
i
=
1
:
model
.
nbStates
%Compute covariance
DataTmp
=
Data

repmat
(
model
.
Mu
(:,
i
),
1
,
nbData
);
S
(:,:,
i
)
=
DataTmp
*
diag
(
GAMMA2
(
i
,:))
*
DataTmp
'
+
eye
(
model
.
nbVar
)
*
diagRegularizationFactor
;
%HDGMM update
, see Eq. (2.2.2) in doc/TechnicalReport.pdf
%HDGMM update
[
V
,
D
]
=
eig
(
S
(:,:,
i
));
[
~
,
id
]
=
sort
(
diag
(
D
),
'descend'
);
% model.D(:,:,i) = D(id(1:model.nbFA), id(1:model.nbFA));
...
...
@@ 40,7 +79,7 @@ for nbIter=1:nbMaxSteps
model
.
D
(:,:,
i
)
=
diag
([
d
(
id
(
1
:
model
.
nbFA
));
repmat
(
mean
(
d
(
id
(
model
.
nbFA
+
1
:
end
))),
model
.
nbVar

model
.
nbFA
,
1
)]);
model
.
V
(:,:,
i
)
=
V
(:,
id
);
%Reconstruct Sigma
, see Eq. (2.2.1) in doc/TechnicalReport.pdf
%Reconstruct Sigma
model
.
Sigma
(:,:,
i
)
=
model
.
V
(:,:,
i
)
*
model
.
D
(:,:,
i
)
*
model
.
V
(:,:,
i
)
'
+
eye
(
model
.
nbVar
)
*
diagRegularizationFactor
;
end
...
...
@@ 57,9 +96,9 @@ end
disp
([
'The maximum number of '
num2str
(
nbMaxSteps
)
' EM iterations has been reached.'
]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function
[
Lik
,
GAMMA
]
=
computeGamma
(
Data
,
model
)
%See Eq. (2.0.5) in doc/TechnicalReport.pdf
Lik
=
zeros
(
model
.
nbStates
,
size
(
Data
,
2
));
for
i
=
1
:
model
.
nbStates
Lik
(
i
,:)
=
model
.
Priors
(
i
)
*
gaussPDF
(
Data
,
model
.
Mu
(:,
i
),
model
.
Sigma
(:,:,
i
));
...
...
m_fcts/EM_MFA.m
View file @
2cb8784d
function
[
model
,
GAMMA2
]
=
EM_MFA
(
Data
,
model
)
%EM for Mixture of factor analysis
%Implementation inspired by "Parsimonious Gaussian Mixture Models" by McNicholas and Murphy, Appendix 8, p.17 (UUU version)
%Sylvain Calinon, 2015
% EM for Mixture of factor analysis (implementation based on "Parsimonious Gaussian
% Mixture Models" by McNicholas and Murphy, Appendix 8, p.17, UUU version).
%
% Writing code takes time. Polishing it and making it available to others takes longer!
% If some parts of the code were useful for your research of for a better understanding
% of the algorithms, please reward the authors by citing the related publications,
% and consider making your own research available in this way.
%
% @article{Calinon15,
% author="Calinon, S.",
% title="A Tutorial on TaskParameterized Movement Learning and Retrieval",
% journal="Intelligent Service Robotics",
% year="2015"
% }
%
% Copyright (c) 2015 Idiap Research Institute, http://idiap.ch/
% Written by Sylvain Calinon, http://calinon.ch/
%
% This file is part of PbDlib, http://www.idiap.ch/software/pbdlib/
%
% PbDlib is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License version 3 as
% published by the Free Software Foundation.
%
% PbDlib is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with PbDlib. If not, see <http://www.gnu.org/licenses/>.
%Parameters of the EM iterations
nbMinSteps
=
5
;
%Minimum number of iterations allowed
...
...
@@ 9,7 +38,7 @@ nbMaxSteps = 100; %Maximum number of iterations allowed
maxDiffLL
=
1E4
;
%Likelihood increase threshold to stop the algorithm
nbData
=
size
(
Data
,
2
);
diagRegularizationFactor
=
1E6
;
%
Regularization term is optional, see Eq. (2.1.2) in doc/TechnicalReport.pdf
diagRegularizationFactor
=
1E6
;
%
Optional regularization term
% %Circular initialization of the MFA parameters
% Itmp = eye(model.nbVar)*1E2;
...
...
@@ 22,7 +51,7 @@ for i=1:model.nbStates
[
V
,
D
]
=
eig
(
model
.
Sigma
(:,:,
i
)

model
.
P
(:,:,
i
));
[
~
,
id
]
=
sort
(
diag
(
D
),
'descend'
);
V
=
V
(:,
id
)
*
D
(
id
,
id
)
.^.
5
;
model
.
L
(:,:,
i
)
=
V
(:,
1
:
model
.
nbFA
);
%>Sigma=LL'+P
model
.
L
(:,:,
i
)
=
V
(:,
1
:
model
.
nbFA
);
%>
Sigma=LL'+P
end
for
nbIter
=
1
:
nbMaxSteps
for
i
=
1
:
model
.
nbStates
...
...
@@ 41,26 +70,26 @@ for nbIter=1:nbMaxSteps
GAMMA2
=
GAMMA
.
/
repmat
(
sum
(
GAMMA
,
2
),
1
,
nbData
);
%Mstep
%Update Priors
, see Eq. (2.2.10) in doc/TechnicalReport.pdf
%Update Priors
model
.
Priors
=
sum
(
GAMMA
,
2
)
/
nbData
;
%Update Mu
, see Eq. (2.2.11) in doc/TechnicalReport.pdf
%Update Mu
model
.
Mu
=
Data
*
GAMMA2
'
;
%Update factor analysers parameters
for
i
=
1
:
model
.
nbStates
%Compute covariance
, see Eq. (2.2.15) in doc/TechnicalReport.pdf
%Compute covariance
DataTmp
=
Data

repmat
(
model
.
Mu
(:,
i
),
1
,
nbData
);
S
(:,:,
i
)
=
DataTmp
*
diag
(
GAMMA2
(
i
,:))
*
DataTmp
'
+
eye
(
model
.
nbVar
)
*
diagRegularizationFactor
;
%Update B
, see Eq. (2.2.16) in doc/TechnicalReport.pdf
%Update B
B
(:,:,
i
)
=
model
.
L
(:,:,
i
)
' / (model.L(:,:,i) * model.L(:,:,i)'
+
model
.
P
(:,:,
i
));
%Update Lambda
, see Eq. (2.2.12) in doc/TechnicalReport.pdf
%Update Lambda
model
.
L
(:,:,
i
)
=
S
(:,:,
i
)
*
B
(:,:,
i
)
' / (eye(model.nbFA)  B(:,:,i) * model.L(:,:,i) + B(:,:,i) * S(:,:,i) * B(:,:,i)'
);
%Update Psi
, see Eq. (2.2.13) in doc/TechnicalReport.pdf
%Update Psi
model
.
P
(:,:,
i
)
=
diag
(
diag
(
S
(:,:,
i
)

model
.
L
(:,:,
i
)
*
B
(:,:,
i
)
*
S
(:,:,
i
)))
+
eye
(
model
.
nbVar
)
*
diagRegularizationFactor
;
%Reconstruct Sigma
, see Eq. (2.2.4) in doc/TechnicalReport.pdf
%Reconstruct Sigma
model
.
Sigma
(:,:,
i
)
=
model
.
L
(:,:,
i
)
*
model
.
L
(:,:,
i
)
'
+
model
.
P
(:,:,
i
);
end
%Compute average loglikelihood
...
...
@@ 76,9 +105,9 @@ end
disp
([
'The maximum number of '
num2str
(
nbMaxSteps
)
' EM iterations has been reached.'
]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function
[
Lik
,
GAMMA
]
=
computeGamma
(
Data
,
model
)
%See Eq. (2.2.9) in doc/TechnicalReport.pdf
Lik
=
zeros
(
model
.
nbStates
,
size
(
Data
,
2
));
for
i
=
1
:
model
.
nbStates
Lik
(
i
,:)
=
model
.
Priors
(
i
)
*
gaussPDF
(
Data
,
model
.
Mu
(:,
i
),
model
.
Sigma
(:,:,
i
));
...
...
m_fcts/EM_MPPCA.m
View file @
2cb8784d
function
[
model
,
GAMMA2
]
=
EM_MPPCA
(
Data
,
model
)
%EM for mixture of probabilistic principal component analyzers,
%inspired by "Mixtures of Probabilistic Principal Component Analysers" by Michael E. Tipping and Christopher M. Bishop
%Sylvain Calinon, 2015
% EM for mixture of probabilistic principal component analyzers (implementation based on
% "Mixtures of Probabilistic Principal Component Analysers" by Michael E. Tipping and Christopher M. Bishop)
%
% Writing code takes time. Polishing it and making it available to others takes longer!
% If some parts of the code were useful for your research of for a better understanding
% of the algorithms, please reward the authors by citing the related publications,
% and consider making your own research available in this way.
%
% @article{Calinon15,
% author="Calinon, S.",
% title="A Tutorial on TaskParameterized Movement Learning and Retrieval",
% journal="Intelligent Service Robotics",
% year="2015"
% }
%
% Copyright (c) 2015 Idiap Research Institute, http://idiap.ch/
% Written by Sylvain Calinon, http://calinon.ch/
%
% This file is part of PbDlib, http://www.idiap.ch/software/pbdlib/
%
% PbDlib is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License version 3 as
% published by the Free Software Foundation.
%
% PbDlib is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with PbDlib. If not, see <http://www.gnu.org/licenses/>.
%Parameters of the EM iterations
nbMinSteps
=
5
;
%Minimum number of iterations allowed
...
...
@@ 29,29 +58,29 @@ for nbIter=1:nbMaxSteps
GAMMA2
=
GAMMA
.
/
repmat
(
sum
(
GAMMA
,
2
),
1
,
nbData
);
%Mstep
%Update Priors
, see Eq. (2.2.10) in doc/TechnicalReport.pdf
%Update Priors
model
.
Priors
=
sum
(
GAMMA
,
2
)
/
nbData
;
%Update Mu
, see Eq. (2.2.11) in doc/TechnicalReport.pdf
%Update Mu
model
.
Mu
=
Data
*
GAMMA2
'
;
%Update factor analyser params
for
i
=
1
:
model
.
nbStates
%Compute covariance
, see Eq. (2.2.19) in doc/TechnicalReport.pdf
%Compute covariance
DataTmp
=
Data

repmat
(
model
.
Mu
(:,
i
),
1
,
nbData
);
S
(:,:,
i
)
=
DataTmp
*
diag
(
GAMMA2
(
i
,:))
*
DataTmp
'
+
eye
(
model
.
nbVar
)
*
diagRegularizationFactor
;
%Update M
, see Eq. (2.2.20) in doc/TechnicalReport.pdf
%Update M
M
=
eye
(
model
.
nbFA
)
*
model
.
o
(
i
)
+
model
.
L
(:,:,
i
)
'
*
model
.
L
(:,:,
i
);
%Update Lambda
, see Eq. (2.2.17) in doc/TechnicalReport.pdf
%Update Lambda
Lnew
=
S
(:,:,
i
)
*
model
.
L
(:,:,
i
)
/
(
eye
(
model
.
nbFA
)
*
model
.
o
(
i
)
+
M
\
model
.
L
(:,:,
i
)
'
*
S
(:,:,
i
)
*
model
.
L
(:,:,
i
));
%Update of sigma^2
, see Eq. (2.2.21) in doc/TechnicalReport.pdf
%Update of sigma^2
model
.
o
(
i
)
=
trace
(
S
(:,:,
i
)

S
(:,:,
i
)
*
model
.
L
(:,:,
i
)
/
M
*
Lnew
'
)
/
model
.
nbVar
;
model
.
L
(:,:,
i
)
=
Lnew
;
%Update Psi
, see Eq. (2.2.18) in doc/TechnicalReport.pdf
%Update Psi
model
.
P
(:,:,
i
)
=
eye
(
model
.
nbVar
)
*
model
.
o
(
i
);
%Reconstruct Sigma
, see Eq. (2.2.4) in doc/TechnicalReport.pdf
%Reconstruct Sigma
model
.
Sigma
(:,:,
i
)
=
model
.
L
(:,:,
i
)
*
model
.
L
(:,:,
i
)
'
+
model
.
P
(:,:,
i
);
end
...
...
@@ 68,9 +97,9 @@ end
disp
([
'The maximum number of '
num2str
(
nbMaxSteps
)
' EM iterations has been reached.'
]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function
[
Lik
,
GAMMA
]
=
computeGamma
(
Data
,
model
)
%See Eq. (2.0.5) in doc/TechnicalReport.pdf
Lik
=
zeros
(
model
.
nbStates
,
size
(
Data
,
2
));
for
i
=
1
:
model
.
nbStates
Lik
(
i
,:)
=
model
.
Priors
(
i
)
*
gaussPDF
(
Data
,
model
.
Mu
(:,
i
),
model
.
Sigma
(:,:,
i
));
...
...
m_fcts/EM_TPGMM.m
View file @
2cb8784d
...
...
@@ 3,21 +3,35 @@ function model = EM_TPGMM(Data, model)
% The approach allows the modulation of the centers and covariance matrices of the Gaussians with respect to
% external parameters represented in the form of candidate coordinate systems.
%
% Author: Sylvain Calinon, 2014
% http://programmingbydemonstration.org/SylvainCalinon
% Writing code takes time. Polishing it and making it available to others takes longer!
% If some parts of the code were useful for your research of for a better understanding
% of the algorithms, please reward the authors by citing the related publications,
% and consider making your own research available in this way.
%
% This source code is given for free! In exchange, I would be grateful if you cite
% the following reference in any academic publication that uses this code or part of it:
%
% @inproceedings{Calinon14ICRA,
% author="Calinon, S. and Bruno, D. and Caldwell, D. G.",
% title="A taskparameterized probabilistic model with minimal intervention control",
% booktitle="Proc. {IEEE} Intl Conf. on Robotics and Automation ({ICRA})",
% year="2014",
% month="MayJune",
% address="Hong Kong, China",
% pages="33393344"
% @article{Calinon15,
% author="Calinon, S.",
% title="A Tutorial on TaskParameterized Movement Learning and Retrieval",
% journal="Intelligent Service Robotics",
% year="2015"
% }
%
% Copyright (c) 2015 Idiap Research Institute, http://idiap.ch/
% Written by Sylvain Calinon, http://calinon.ch/
%
% This file is part of PbDlib, http://www.idiap.ch/software/pbdlib/
%
% PbDlib is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License version 3 as
% published by the Free Software Foundation.
%
% PbDlib is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with PbDlib. If not, see <http://www.gnu.org/licenses/>.
%Parameters of the EM algorithm
nbMinSteps
=
5
;
%Minimum number of iterations allowed
...
...
@@ 25,8 +39,8 @@ nbMaxSteps = 100; %Maximum number of iterations allowed
maxDiffLL
=
1E5
;
%Likelihood increase threshold to stop the algorithm
nbData
=
size
(
Data
,
3
);
diagRegularizationFactor
=
1E5
;
%
Regularization term is optional, see Eq. (2.1.2) in doc/TechnicalReport.pdf
%diagRegularizationFactor = 0; %
Regularization term is optional, see Eq. (2.1.2) in doc/TechnicalReport.pdf
diagRegularizationFactor
=
1E5
;
%
Optional regularization term
%diagRegularizationFactor = 0; %
Optional regularization term
for
nbIter
=
1
:
nbMaxSteps
fprintf
(
'.'
);
...
...
@@ 38,7 +52,7 @@ for nbIter=1:nbMaxSteps
%Mstep
for
i
=
1
:
model
.
nbStates
%Update Priors
, see Eq. (6.0.2) in doc/TechnicalReport.pdf
%Update Priors
model
.
Priors
(
i
)
=
sum
(
sum
(
GAMMA
(
i
,:)))
/
nbData
;
for
m
=
1
:
model
.
nbFrames
...
...
@@ 46,10 +60,10 @@ for nbIter=1:nbMaxSteps
DataMat
=
[];
DataMat
(
1
:
model
.
nbVars
(
m
),:)
=
Data
(
1
:
model
.
nbVars
(
m
),
m
,:);
%Update Mu
, see Eq. (6.0.3) in doc/TechnicalReport.pdf
%Update Mu
model
.
Mu
(
1
:
model
.
nbVars
(
m
),
m
,
i
)
=
DataMat
*
GAMMA2
(
i
,:)
'
;
%Update Sigma (regularization term is optional)
, see Eq. (6.0.4) in doc/TechnicalReport.pdf
%Update Sigma (regularization term is optional)
DataTmp
=
DataMat

repmat
(
model
.
Mu
(
1
:
model
.
nbVars
(
m
),
m
,
i
),
1
,
nbData
);
model
.
Sigma
(
1
:
model
.
nbVars
(
m
),
1
:
model
.
nbVars
(
m
),
m
,
i
)
=
DataTmp
*
diag
(
GAMMA2
(
i
,:))
*
DataTmp
'
+
eye
(
model
.
nbVars
(
m
))
*
diagRegularizationFactor
;
end
...
...
@@ 68,9 +82,9 @@ end
disp
([
'The maximum number of '
num2str
(
nbMaxSteps
)
' EM iterations has been reached.'
]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function
[
Lik
,
GAMMA
,
GAMMA0
]
=
computeGamma
(
Data
,
model
)
%See Eq. (6.0.1) in doc/TechnicalReport.pdf
nbData
=
size
(
Data
,
3
);
Lik
=
ones
(
model
.
nbStates
,
nbData
);
GAMMA0
=
zeros
(
model
.
nbStates
,
model
.
nbFrames
,
nbData
);
...
...
@@ 85,7 +99,3 @@ for i=1:model.nbStates
end
GAMMA
=
Lik
.
/
repmat
(
sum
(
Lik
,
1
)
+
realmin
,
size
(
Lik
,
1
),
1
);
end
m_fcts/EM_stdPGMM.m
View file @
2cb8784d
function
[
model
,
s
,
LL
]
=
EM_stdPGMM
(
s
,
model
)
% Training of a parametric Gaussian mixture model (PGMM) with expectationmaximization (EM).
% The implementation follows the approach described by Wilson and Bobick (1999) "Parametric Hidden Markov
% Models for Gesture Recognition", IEEE Trans. on Pattern Analysis and Machine Intelligence, with an
% implementation applied to the special case of Gaussian mixture models (GMM).
% The approach is inspired by Wilson and Bobick (1999), with an implementation applied to
% the special case of Gaussian mixture models (GMM).
%
% Author: Sylvain Calinon, 2013
% http://programmingbydemonstration.org/SylvainCalinon
% Writing code takes time. Polishing it and making it available to others takes longer!
% If some parts of the code were useful for your research of for a better understanding
% of the algorithms, please reward the authors by citing the related publications,
% and consider making your own research available in this way.
%
% This source code is given for free! In exchange, I would be grateful if you cite
% the following references in any academic publication that uses this code or part of it:
%
% @inproceedings{Calinon12Hum,
% author="Calinon, S. and Li, Z. and Alizadeh, T. and Tsagarakis, N. G. and Caldwell, D. G.",
% title="Statistical dynamical systems for skills acquisition in humanoids",
% booktitle="Proc. {IEEE} Intl Conf. on Humanoid Robots ({H}umanoids)",
% year="2012",
% address="Osaka, Japan"
% @article{Calinon15,
% author="Calinon, S.",
% title="A Tutorial on TaskParameterized Movement Learning and Retrieval",
% journal="Intelligent Service Robotics",
% year="2015"
% }
%
% @article{Wilson99,
% author="Wilson, A. D. and Bobick, A. F.",
% title="Parametric Hidden {M}arkov Models for Gesture Recognition",
...
...
@@ 27,8 +25,23 @@ function [model, s, LL] = EM_stdPGMM(s, model)
% pages="884900"
% }
%
% The first reference presents an implementation of the approach described in the second reference, and
% applies it to the special case of Gaussian mixture model (GMM).
% Copyright (c) 2015 Idiap Research Institute, http://idiap.ch/
% Written by Sylvain Calinon, http://calinon.ch/
%
% This file is part of PbDlib, http://www.idiap.ch/software/pbdlib/
%
% PbDlib is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License version 3 as
% published by the Free Software Foundation.
%
% PbDlib is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with PbDlib. If not, see <http://www.gnu.org/licenses/>.
%Parameters of the EM algorithm
nbMinSteps
=
10
;
%Minimum number of iterations allowed
...
...
@@ 54,19 +67,19 @@ for nbIter=1:nbMaxSteps
%MSTEP
for
i
=
1
:
model
.
nbStates
%Update Priors
, see Eq. (7.1.5) in doc/TechnicalReport.pdf
%Update Priors
model
.
Priors
(
i
)
=
sum
(
GAMMA
(
i
,:))/
nbData
;
%Update Zmu
, see Eq. (7.1.6) in doc/TechnicalReport.pdf
%Update Zmu
model
.
ZMu
(:,:,
i
)
=
zeros
(
model
.
nbVar
,
nbVarParams
);
sumTmp
=
zeros
(
nbVarParams
,
nbVarParams
);
for
n
=
1
:
nbSamples
model
.
ZMu
(:,:,
i
)
=
model
.
ZMu
(:,:,
i
)
+
(
s
(
n
)
.
Data
*
diag
(
s
(
n
)
.
GAMMA
(
i
,:))
*
repmat
(
s
(
n
)
.
OmegaMu
'
,
s
(
n
)
.
nbData
,
1
));
sumTmp
=
sumTmp
+
(
s
(
n
)
.
OmegaMu
*
s
(
n
)
.
OmegaMu
'
)
*
(
sum
(
s
(
n
)
.
GAMMA
(
i
,:)));
end
model
.
ZMu
(:,:,
i
)
=
model
.
ZMu
(:,:,
i
)
*
pinv
(
sumTmp
+
eye
(
nbVarParams
)
*
diagRegularizationFactor
);
%Eq. (6) Wilson and Bobick
model
.
ZMu
(:,:,
i
)
=
model
.
ZMu
(:,:,
i
)
*
pinv
(
sumTmp
+
eye
(
nbVarParams
)
*
diagRegularizationFactor
);
%Update Sigma
, see Eq. (7.1.7) in doc/TechnicalReport.pdf
%Update Sigma
model
.
Sigma
(:,:,
i
)
=
zeros
(
model
.
nbVar
);
sumTmp
=
0
;
for
n
=
1
:
nbSamples
...
...
@@ 95,7 +108,6 @@ for nbIter=1:nbMaxSteps
end
end
disp
([
'The maximum number of '
num2str
(
nbMaxSteps
)
' EM iterations has been reached.'
]);
end
...
...
@@ 121,7 +133,3 @@ for n=1:nbSamples
nTmp
=
nTmp
+
size
(
s
(
n
)
.
GAMMA
,
2
);
end
end
m_fcts/EM_tensorGMM.m
View file @
2cb8784d
...
...
@@ 3,21 +3,35 @@ function [model, GAMMA0, GAMMA2] = EM_tensorGMM(Data, model)
% The approach allows the modulation of the centers and covariance matrices of the Gaussians with respect to
% external parameters represented in the form of candidate coordinate systems.
%
% Author: Sylvain Calinon, 2014
% http://programmingbydemonstration.org/SylvainCalinon
% Writing code takes time. Polishing it and making it available to others takes longer!
% If some parts of the code were useful for your research of for a better understanding
% of the algorithms, please reward the authors by citing the related publications,
% and consider making your own research available in this way.
%
% This source code is given for free! In exchange, I would be grateful if you cite
% the following reference in any academic publication that uses this code or part of it:
%
% @inproceedings{Calinon14ICRA,
% author="Calinon, S. and Bruno, D. and Caldwell, D. G.",
% title="A taskparameterized probabilistic model with minimal intervention control",
% booktitle="Proc. {IEEE} Intl Conf. on Robotics and Automation ({ICRA})",
% year="2014",
% month="MayJune",
% address="Hong Kong, China",
% pages="33393344"
% @article{Calinon15,
% author="Calinon, S.",
% title="A Tutorial on TaskParameterized Movement Learning and Retrieval",
% journal="Intelligent Service Robotics",
% year="2015"
% }
%
% Copyright (c) 2015 Idiap Research Institute, http://idiap.ch/
% Written by Sylvain Calinon, http://calinon.ch/
%
% This file is part of PbDlib, http://www.idiap.ch/software/pbdlib/
%
% PbDlib is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License version 3 as
% published by the Free Software Foundation.
%
% PbDlib is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with PbDlib. If not, see <http://www.gnu.org/licenses/>.
%Parameters of the EM algorithm
nbMinSteps
=
5
;
%Minimum number of iterations allowed
...
...
@@ 25,8 +39,8 @@ nbMaxSteps = 100; %Maximum number of iterations allowed
maxDiffLL
=
1E5
;
%Likelihood increase threshold to stop the algorithm
nbData
=
size
(
Data
,
3
);