1

A Bayesian Approach to Policy Recognition
and State Representation Learning

arXiv:1605.01278v4 [stat.ML] 4 Aug 2017

Adrian Šošić, Abdelhak M. Zoubir and Heinz Koeppl
Abstract—Learning from demonstration (LfD) is the process of building behavioral models of a task from demonstrations provided by
an expert. These models can be used e.g. for system control by generalizing the expert demonstrations to previously unencountered
situations. Most LfD methods, however, make strong assumptions about the expert behavior, e.g. they assume the existence of a
deterministic optimal ground truth policy or require direct monitoring of the expert’s controls, which limits their practical use as part of a
general system identification framework. In this work, we consider the LfD problem in a more general setting where we allow for
arbitrary stochastic expert policies, without reasoning about the optimality of the demonstrations. Following a Bayesian methodology,
we model the full posterior distribution of possible expert controllers that explain the provided demonstration data. Moreover, we show
that our methodology can be applied in a nonparametric context to infer the complexity of the state representation used by the expert,
and to learn task-appropriate partitionings of the system state space.
Index Terms—learning from demonstration, policy recognition, imitation learning, Bayesian nonparametric modeling, Markov chain
Monte Carlo, Gibbs sampling, distance dependent Chinese restaurant process

F

1

L

I NTRODUCTION

(LfD) has become a
viable alternative to classical reinforcement learning as
a new data-driven learning paradigm for building behavioral models based on demonstration data. By exploiting
the domain knowledge provided by an expert demonstrator,
LfD-built models can focus on the relevant parts of a system’s state space [1] and hence avoid the need of tedious exploration steps performed by reinforcement learners, which
often require an impractically high number of interactions
with the system environment [2] and always come with the
risk of letting the system run into undesired or unsafe states
[3]. In addition to that, LfD-built models have been shown
to outperform the expert in several experiments [4], [5], [6].
However, most existing LfD methods come with strong
requirements that limit their practical use in real-world scenarios. In particular, they often require direct monitoring of
the expert’s controls (e.g. [5], [7], [8]) which is possible only
under laboratory-like conditions, or they need to interact
with the target system via a simulator, if not by controlling
the system directly (e.g. [9]). Moreover, many methods are
restricted to problems with finite state spaces (e.g. [10]), or
they compute only point estimates of the relevant system
parameters without providing any information about their
level of confidence (e.g. [9], [11], [12]). Last but not least, the
expert is typically assumed to follow an optimal deterministic policy (e.g. [13]) or to at least approximate one, based
on some presupposed degree of confidence in the optimal
•
•
•

EARNING FROM DEMONSTRATION

Adrian Šošić is a member of the Signal Processing Group and an associate
member of the Bioinspired Communication Systems Lab, Technische Universität Darmstadt, Germany. E-mail: adrian.sosic@spg.tu-darmstadt.de
Abdelhak M. Zoubir is the head of the Signal Processing Group, Technische
Universität Darmstadt, Germany. E-mail: zoubir@spg.tu-darmstadt.de
Heinz Koeppl is the head of the Bioinspired Communication Systems Lab
and a member of the Centre for Cognitive Science, Technische Universität
Darmstadt, Germany. E-mail: heinz.koeppl@bcs.tu-darmstadt.de

behavior (e.g. [14]). While such an assumption may be
reasonable in some situations (e.g. for problems in robotics
involving a human demonstrator [1]), it is not appropriate
in many others, such as in multi-agent environments, where
an optimal deterministic policy often does not exist [15].
In fact, there are many situations in which the assumption
of a deterministic expert behavior is violated. In a more
general system identification setting, our goal could be, for
instance, to detect the deviation of an agent’s policy from
its known nominal behavior, e.g. for the purpose of fault or
fraud detection (note that the term “expert” is slightly misleading in this context). Also, there are situations in which
we might not want to reason about the optimality of the
demonstrations; for instance, when studying the exploration
strategy of an agent who tries to model its environment (or
the reactions of other agents [16]) by randomly triggering
different events. In all these cases, existing LfD methods
can at best approximate the behavior of the expert as they
presuppose the existence of some underlying deterministic
ground truth policy.
In this work, we present a novel approach to LfD in order
to address the above-mentioned shortcomings of existing
methods. Central to our work is the problem of policy
recognition, that is, extracting the (possibly stochastic and
non-optimal) policy of a system from observations of its
behavior. Taking a general system identification view on the
problem, our goal is herein to make as few assumptions
about the expert behavior as possible. In particular, we consider the whole class of stochastic expert policies, without
ever reasoning about the optimality of the demonstrations.
As a result of this, our hypothesis space is not restricted to a
certain class of ground truth policies, such as deterministic
or softmax policies (c.f. [14]). This is in contrast to inverse
reinforcement learning approaches (see Section 1.2), which
interpret the observed demonstrations as the result of some
preceding planing procedure conducted by the expert which

2

they try to invert. In the above-mentioned case of fault
detection, for example, such an inversion attempt will generally fail since the demonstrated behavior can be arbitrarily
far from optimal, which renders an explanation of the data
in terms of a simple reward function impossible.
Another advantage of our problem formulation is that
the resulting inference machinery is entirely passive, in the
sense that we require no active control of the target system
nor access to the action sequence performed by the expert.
Accordingly, our method is applicable to a broader range of
problems than targeted by most existing LfD frameworks
and can be used for system identification in cases where
we cannot interact with the target system. However, our
objective in this paper is twofold: we not only attempt to
answer the question how the expert performs a given task
but also to infer which information is used by the expert to
solve it. This knowledge is captured in the form of a joint
posterior distribution over possible expert state representations and corresponding state controllers. As the complexity
of the expert’s state representation is unknown a priori,
we finally present a Bayesian nonparametric approach to
explore the underlying structure of the system space based
on the available demonstration data.
1.1 Problem statement
Given a set of expert demonstrations in the form of a system
trajectory s = (s1 , s2 , . . . , sT ) ∈ S T of length T , where
S denotes the system state space, our goal is to determine
the latent control policy used by the expert to generate the
state sequence.1 We formalize this problem as a discretetime decision-making process (i.e. we assume that the expert
executes exactly one control action per trajectory state) and
adopt the Markov decision process (MDP) formalism [17]
as the underlying framework describing the dynamics of
our system. More specifically, we consider a reduced MDP
(S, A, T , π) which consists of a countable or uncountable
system state space S , a finite set of actions A containing
|A| elements, a transition model T : S × S × A → R≥0
where T (s0 | s, a) denotes the probability (density) assigned
to the event of reaching state s0 after taking action a in
state s, and a policy π modeling the expert’s choice of
actions.2 In the following, we assume that the expert policy
is parametrized by a parameter ω ∈ Ω, which we call
the global control parameter of the system, and we write
π(a | s, ω), π : A × S × Ω → [0, 1], to denote the expert’s
local policy (i.e. the distribution of actions a played by the
expert) at any given state s under ω . The set Ω is called the
parameter space of the policy, which specifies the class of
feasible action distributions. The specific form of Ω will be
discussed later.
Using a parametric description for π is convenient as
it shifts the recognition task from determining the possibly
infinite set of local policies at all states in S to inferring the
posterior distribution p(ω | s), which contains all information that is relevant for predicting the expert behavior,
Z
p(a | s∗ , s) =
π(a | s∗ , ω)p(ω | s) dω.
Ω

1. The generalization to multiple trajectories is straightforward as
they are conditionally independent given the system parameters.
2. This reduced model is sometimes referred to as an MDP\R (see
e.g. [9], [18], [19]) to emphasize the nonexistence of a reward function.

Herein, s∗ ∈ S is some arbitrary query point and p(a | s∗ , s)
is the corresponding predictive action distribution. Since the
local policies are coupled through the global control parameter ω as indicated by the above integral equation, inferring
ω means not only to determine the individual local policies
but also their spatial dependencies. Consequently, learning
the structure of ω from demonstration data can be also
interpreted as learning a suitable state representation for
the task performed by the expert. This relationship will be
discussed in detail in the forthcoming sections. In Section 3,
we further extend this reasoning to a nonparametric policy
model whose hypothesis class finally covers all stochastic
policies on S .
For the remainder of this paper, we make the common
assumptions that the transition model T as well as the
system state space S and the action set A are known. The
assumption of knowing S follows naturally because we
already assumed that we can observe the expert acting in S .
In the proposed Bayesian framework, the latter assumption
can be easily relaxed by considering noisy or incomplete trajectory data. However, as this would not provide additional
insights into the main principles of our method, we do not
consider such an extension in this work.
The assumption of knowing the transition dynamics T is
a simplifying one but prevents us from running into model
identifiability problems: if we do not constrain our system
transition model in some reasonable way, any observed state
transition in S could be trivially explained by a corresponding local adaptation of the assumed transition model T
and, thus, there would be little hope to extract the true expert policy from the demonstration data. Assuming a fixed
transition model is the easiest way to resolve this model
ambiguity. However, there are alternatives which we leave
for future work, for example, using a parametrized family of
transition models for joint inference. This extension can be
integrated seamlessly into our Bayesian framework and is
useful in cases where we can constrain the system dynamics
in a natural way, e.g. when modeling physical processes.
Also, it should be mentioned that we can tolerate deviations
from the true system dynamics as long as our model T is
sufficiently accurate to extract information about the expert
action sequence locally, because our inference algorithm
naturally processes the demonstration data piece-wise in
the form of one-step state transitions {(st , st+1 )} (see algorithmic details in Section 2 and results in Section 4.2). This
is in contrast to planning-based approaches, where small
modeling errors in the dynamics can accumulate and yield
consistently wrong policy estimates [8], [20].
The requirement of knowing the action set A is less
stringent: if A is unknown a priori, we can still assume a potentially rich class of actions, as long as the transition model
can provide the corresponding dynamics (see example in
Section 4.2). For instance, we might be able to provide a
model which describes the movement of a robotic arm even
if the maximum torque that can be generated by the system
is unknown. Figuring out which of the hypothetical actions
are actually performed by the expert and, more importantly,
how they are used in a given context, shall be the task of
our inference algorithm.

3

1.2

Related work

The idea of learning from demonstration has now been
around for several decades. Most of the work on LfD has
been presented by the robotics community (see [1] for a
survey), but recent advances in the field have triggered developments in other research areas, such as cognitive science
[21] and human-machine interaction [22]. Depending on the
setup, the problem is referred to as imitation learning [23],
apprenticeship learning [9], inverse reinforcement learning
[13], inverse optimal control [24], preference elicitation [21],
plan recognition [25] or behavioral cloning [5]. Most LfD
models can be categorized as intentional models (with
inverse reinforcement learning models as the primary example), or sub-intentional models (e.g. behavioral cloning
models). While the latter class only predicts an agent’s
behavior via a learned policy representation, intentional
models (additionally) attempt to capture the agent’s beliefs
and intentions, e.g. in the form of a reward function. For this
reason, intentional models are often reputed to have better
generalization abilities3 ; however, they typically require a
certain amount of task-specific prior knowledge in order
to resolve the ambiguous relationship between intention
and behavior, since there are often many ways to solve
a certain task [13]. Also, albeit being interesting from a
psychological point of view [21], intentional models target
a much harder problem than what is actually required in
many LfD scenarios. For instance, it is not necessary to
understand an agent’s intention if we only wish to analyze
its behavior locally.
Answering the question whether or not an intentionbased modeling of the LfD problem is advantageous, is out
of the scope of this paper; however, we point to the comprehensive discussion in [26]. Rather, we present a hybrid
solution containing both intentional and sub-intentional elements. More specifically, our method does not explicitly capture the expert’s goals in the form of a reward function but
infers a policy model directly from the demonstration data;
nonetheless, the presented algorithm learns a task-specific
representation of the system state space which encodes the
structure of the underlying control problem to facilitate the
policy prediction task. An early version of this idea can
be found in [27], where the authors proposed a simple
method to partition a system’s state space into a set of socalled control situations to learn a global system controller
based on a small set of informative states. However, their
framework does not incorporate any demonstration data
and the proposed partitioning is based on heuristics. A more
sophisticated partitioning approach utilizing expert demonstrations is shown in [11]; yet, the proposed expectationmaximization framework applies to deterministic policies
and finite state spaces only.
The closest methods to ours can be probably found in
[19] and [10]. The authors of [19] presented a nonparametric
inverse reinforcement learning approach to cluster the expert data based on a set of learned subgoals encoded in the
form of local rewards. Unfortunately, the required subgoal
3. The rationale behind this is that an agent’s intention is always
specific to the task being performed and can hence serve as a compact
description of it [13]. However, if the intention of the agent is misunderstood, then also the subsequent generalization step will trivially fail.

assignments are learned only for the demonstration set and,
thus, the algorithm cannot be used for action prediction
at unvisited states unless it is extended with a non-trivial
post-processing step which solves the subgoal assignment
problem. Moreover, the algorithm requires an MDP solver,
which causes difficulties for systems with uncountable state
spaces. The sub-intentional model in [10], on the other hand,
can be used to learn a class of finite state controllers directly
from the expert demonstrations. Like our framework, the
algorithm can handle various kinds of uncertainty about
the data but, again, the proposed approach is limited to
discrete settings. In the context of reinforcement learning,
we further point to the work presented in [28] whose
authors follow a nonparametric strategy similar to ours, to
learn a distribution over predictive state representations for
decision-making.
1.3

Paper outline

The outline of the paper is as follows: In Section 2, we
introduce our parametric policy recognition framework and
derive inference algorithms for both countable and uncountable state spaces. In Section 3, we consider the policy recognition problem from a nonparametric viewpoint and provide insights into the state representation learning problem.
Simulation results are presented in Section 4 and we give
a conclusion of our work in Section 5. In the supplement,
we provide additional simulation results, a note on the
computational complexity of our model, as well as an indepth discussion on the issue of marginal invariance and
the problem of policy prediction in large states spaces.

2
2.1

PARAMETRIC P OLICY R ECOGNITION
Finite state spaces: the static model

First, let us assume that the expert system can be modeled
on a finite state space S and let |S| denote its cardinality.
For notational convenience, we represent both states and
actions by integer values. Starting with the most general
case, we assume that the expert executes an individual
control strategy at each possible system state. Accordingly,
we introduce a set of local control parameters or local con|S|
trollers {θ i }i=1 by which we describe the expert’s choice of
actions. More specifically, we model the executed actions
as categorical random variables and let the j th element of
θ i represent the probability that the expert chooses action
j at state i. Consequently, θ i lies in the (|A| − 1)-simplex,
which we denote by the symbol ∆ for brevity of notation, i.e.
θ i ∈ ∆ ⊆ R|A| . Summarizing all local control parameters in
a single matrix, Θ ∈ Ω ⊆ ∆|S| , we obtain the global control
parameter of the system as already introduced in Section 1.1,
which compactly captures the expert behavior. Note that
we denote the global control parameter here by Θ instead
of ω , for reasons that will become clear later. Each action a
is thus characterized by the local policy that is induced by
the control parameter of the underlying state,

π(a | s = i, Θ) = C AT(a | θ i ).
For simplicity, we will write π(a | θ i ) since the state information is used only to select the appropriate local controller.

4

Considering a finite set of actions, it is convenient to
place a symmetric Dirichlet prior on the local control parameters,
pθ (θ i | α) = D IR(θ i | α · 1|A| ),
which forms the conjugate distribution to the categorical
distribution over actions. Here, 1|A| denotes the vector of
all ones of length |A|. The prior is itself parametrized by a
concentration parameter α which can be further described
by a hyperprior pα (α), giving rise to a Bayesian hierarchical
model. For simplicity, we assume that the value of α is fixed
for the remainder of this paper, but the extension to a full
Bayesian treatment is straightforward. The joint distribution
of all remaining model variables is, therefore, given as

p(s, a, Θ | α) = p1 (s1 )
... ×

TY
−1
t=1

|S|
Y

i=1

pθ (θ i | α) . . .

(1)

T (st+1 | st , at )π(at | θ st ),

where a = (a1 , a2 , . . . , aT −1 ) denotes the latent action
sequence taken by the expert and p1 (s1 ) is the initial state
distribution of the system. Throughout the rest of the paper,
we refer to this model as the static model. The corresponding
graphical visualization is depicted in Fig. 1.
2.1.1

Gibbs sampling

Following a Bayesian methodology, our goal is to determine
the posterior distribution p(Θ | s, α), which contains all
information necessary to make predictions about the expert behavior. For the static model in Eq. (1), the required
marginalization of the latent action sequence a can be
computed efficiently because the joint distribution factorizes
over time instants. For the extended models presented in
later sections, however, a direct marginalization becomes
computationally intractable due to the exponential growth
of latent variable configurations. As a solution to this problem, we follow a sampling-based inference strategy which
is later on generalized to more complex settings.
For the simple model described above, we first approximate the joint posterior distribution p(Θ, a | s, α) over both
controllers and actions using a finite number of Q samples,
and then marginalize over a in a second step,
X
p(Θ | s, α) =
p(Θ, a | s, α)
(2)

≈

X
a



st−1

st

st+1

at−1

at

at+1

α

z

Θ

Fig. 1: Graphical model of the policy recognition framework. The underlying dynamical structure is that of an
MDP whose global control parameter Θ is treated as a
random variable with prior distribution parametrized by α.
The indicator node z is used for the clustering model in
Section 2.2. Observed variables are highlighted in gray.
Herein, a−t and Θ−i refer to all actions/controllers except
at and θ i , respectively. The latter of the two expressions
reveals that, in order to sample θ i , we need to consider
only those actions played at the corresponding state i.
Furthermore, the first expression shows that, given Θ, all
actions {at } can be sampled independently of each other.
Therefore, inference can be done in parallel for all θ i . This
can be also seen from the posterior distribution of the global
control parameter, which factorizes over states,

p(Θ | s, a, α) ∝

|S|
Y

i=1

pθ (θ i | α)

Y

t:st =i

π(at | θ i ).

(3)

From the conjugacy of pθ (θ i | α) and π(at | θ i ), it follows
that the posterior over θ i is again Dirichlet distributed with
updated concentration parameter. In particular, denoting by
φi,j the number of times that action j is played at state i for
the current assignment of actions a,
X
φi,j :=
1(at = j),
(4)
t:st =i

and by collecting these quantities in the form of vectors, i.e.
φi := [φi,1 , . . . , φi,|A| ], we can rewrite Eq. (3) as

p(Θ | s, a, α) =

|S|
Y

i=1

D IR(θ i | φi + α · 1|A| ).

(5)

a


Q
Q
X
1
1 X


δ {q} {q} (Θ, a) =
δ {q} (Θ),
Q q=1 Θ a
Q q=1 Θ

where (Θ{q} , a{q} ) ∼ p(Θ, a | s, α), and δx (·) denotes
Dirac’s delta function centered at x. This two-step approach
gives rise to a simple inference procedure since the joint
Q
samples {(Θ{q} , a{q} )}q=1 can be easily obtained from a
Gibbs sampling scheme, i.e. by sampling iteratively from
the following two conditional distributions,

p(at | a−t , s, Θ, α) ∝ T (st+1 | st , at )π(at | θ st ),
Y
p(θ i | Θ−i , s, a, α) ∝ pθ (θ i | α)
π(at | θ i ).
t:st =i

2.1.2 Collapsed Gibbs sampling
Choosing a Dirichlet distribution as prior model for the
local controllers is convenient as it allows us to arrive at
analytic expressions for the conditional distributions that
are required to run the Gibbs sampler. As an alternative,
we can exploit the conjugacy property of pθ (θ i | α) and
π(at | θ i ) to marginalize out the control parameters during
the sampling process, giving rise to a collapsed sampling
scheme. Collapsed sampling is advantageous in two different respects: first, it reduces the total number of variables
to be sampled and, hence, the number of computations
required per Gibbs iteration; second, it increases the mixing
speed of the underlying Markov chain that governs the
sampling process, reducing the correlation of the obtained

5

samples and, with it, the variance of the resulting policy
estimate.
Formally, collapsing means that we no longer approximate the joint distribution p(Θ, a | s, α) as done in Eq. (2),
but instead sample from the marginal density p(a | s, α),
X
p(Θ | s, α) =
p(Θ | s, a, α)p(a | s, α)
a

≈

X
a



1
p(Θ | s, a, α) 
Q
Q

Q
X

q=1

1 X
=
p(Θ | s, a{q} , α),
Q q=1



δa{q} (a)
(6)

where a{q} ∼ p(a | s, α). In contrast to the previous
approach, the target distribution is no longer represented
by a sum of Dirac measures but described by a product of
Dirichlet mixtures (compare Eq. (5)). The required samples
{a{q} } can be obtained from a collapsed Gibbs sampler with
Z
p(s, a, Θ | α) dΘ
p(at | a−t , s, α) ∝
|S|
∆
Z
Y
∝ T (st+1 | st , at )
pθ (θ st | α)
π(at0 | θ st ) dθ st .
∆

t0 :st0 =st

It turns out that the above distribution provides an easy
sampling mechanism since the integral part, when viewed
as a function of action at only, can be identified as the conditional of a Dirichlet-multinomial distribution. This distribution is then reweighted by the likelihood T (st+1 | st , at )
of the observed transition. The final (unnormalized) weights
of the resulting categorical distribution are hence given as

p(at = j | a−t , s, α) ∝ T (st+1 | st , at = j) · (ϕt,j + α), (7)
where ϕt,j counts the number of occurrences of action j
among all actions in a−t played at the same state as at (that
is, st ). Explicitly,
X
ϕt,j :=
1(at0 = j).
t0 :st0 =st
t0 6=t

Note that these values can be also expressed in terms of the
sufficient statistics introduced in the last section,

ϕt,j = φst ,j − 1(at = j).
As before, actions played at different states may be sampled
independently of each other because they are generated by
different local controllers. Consequently, inference about Θ
again decouples for all states.
2.2

Towards large state spaces: a clustering approach

While the methodology introduced so far provides a means
to solve the policy recognition problem in finite state spaces,
the presented approaches quickly become infeasible for
large spaces as, in the continuous limit, the number of
parameters to be learned (i.e. the size of Θ) will grow
unbounded. In that sense, the presented methodology is
prone to overfitting because, for larger problems, we will
never have enough demonstration data to sufficiently cover
the whole system state space. In particular, the static model

S

θ1

θ2
θ3

θ4
θ5

Fig. 2: Schematic illustration of the clustering model. The
state space S is partitioned into a set of clusters {Ck }, each
governed by its own local control parameter θ k .
makes no assumptions about the structure of Θ but treats
all local policies separately (see Eq. (5)); hence, we are not
able to generalize the demonstrated behavior to regions of
the state space that are not directly visited by the expert.
Yet, we would certainly like to predict the expert behavior
also at states for which there is no trajectory data available. Moreover, we should expect a well-designed model to
produce increasingly accurate predictions at regions closer
to the observed trajectories (with the precise definition of
“closeness” being left open for the moment).
A simple way to counteract the overfitting problem, in
general, is to restrict the complexity of a model by limiting
the number of its free parameters. In our case, we can avoid
the parameter space to grow unbounded by considering
only a finite number of local policies that need to be shared
between the states. The underlying assumption is that, at
each state, the expert selects an action according to one
of K local policies, with corresponding control parameters {θ k }K
k=1 . Accordingly, we introduce a set of indicator
|S|
or cluster assignment variables, {zi }i=1 , zi ∈ {1, . . . , K},
which map the states to their local controllers (Fig. 1).
Obviously, such an assignment implies a partitioning of the
state space (Fig. 2), resulting in the following K clusters,

Ck := {i : zi = k},

k ∈ {1, . . . , K}.

Although we motivated the clustering of states by the
problem of overfitting, partitioning a system’s space is not
only convenient from a statistical point of view; mapping
the inference problem down to a lower-dimensional space
is also reasonable for practical reasons as we are typically
interested in understanding an agent’s behavior on a certain
task-appropriate scale. The following paragraphs discuss
these reasons in detail:

• In practice, the observed trajectory data will always be
noisy since we can take our measurements only up to a
certain finite precision. Even though we do not explicitly
consider observation noise in this paper, clustering the data
appears reasonable in order to robustify the model against
small perturbations in our observations.
• Considering the LfD problem from a control perspective,
the complexity of subsequent planning steps can be potentially reduced if the system dynamics can be approximately
described on a lower-dimensional manifold of the state
space, meaning that the system behavior can be well represented by a smaller set of informative states (c.f. finite state

6

controllers [29], control situations [27]). The LfD problem
can then be interpreted as the problem of learning a (nearoptimal) controller based on a small set of local policies that
together provide a good approximation of the global agent
behavior. What remains is the question how we can find
such a representation. The clustering approach described
above offers one possible solution to this problem.

• Finally, in any real setup, it is reasonable to assume that
the expert itself can only execute a finite-precision policy
due to its own limited sensing abilities of the system state
space. Consequently, the demonstrated behavior is going
to be optimal only up to a certain finite precision because
the agent is generally not able to discriminate between
arbitrary small differences of states. An interesting question
in this context is whether we can infer the underlying state
representation of the expert by observing its reactions to the
environment in the form of the resulting state trajectory. We
will discuss this issue in detail in Section 3.
By introducing the cluster assignment variables {zi }, the
joint distribution in Eq. (1) changes into

p(s, a, z, Θ | α) = p1 (s1 )
... ×

TY
−1
t=1

K
Y

k=1

(8)

pθ (θ k | α) . . .

T (st+1 | st , at )π(at | θ zst )pz (z),

where z = (z1 , z2 , . . . , z|S| ) denotes the collection of all
indicator variables and pz (z) is the corresponding prior
distribution to be further discussed in Section 2.2.3. Note
that the static model can be recovered as a special case of
the above when each state describes its own cluster, i.e. by
setting K = |S| and fixing zi = i (hence the name static).
In contrast to the static model, we now require both the
indicator zi and the corresponding control parameter θ zi
in order to characterize the expert’s behavior at a given
state i. Accordingly, the global control parameter of the
model is ω = (Θ, z) with underlying parameter space
Ω ⊆ ∆K ×{1, . . . , K}|S| (see Section 1.1), and our target distribution becomes p(Θ, z | s, α). In what follows, we derive
the Gibbs and the collapsed Gibbs sampler as mechanisms
for approximate inference in this setting.
2.2.1 Gibbs sampling
As shown by the following equations, the expressions for
the conditional distributions over actions and controllers
take a similar form to those of the static model. Here, the
only difference is that we no longer group the actions by
their states but according to their generating local policies
or, equivalently, the clusters {Ck },

p(at | a−t , s, z, Θ, α) ∝ T (st+1 | st , at ) · π(at | θ zst ),
Y
p(θ k | Θ−k , s, a, z, α) ∝ pθ (θ k | α)
π(at | θ k )
t:zst =k

= pθ (θ k | α)

Y

t:st ∈Ck

π(at | θ k ).

The latter expression again takes the form of a Dirichlet
distribution with updated concentration parameter,

p(θ k | Θ−k , s, a, z, α) = D IR(θ k | ξ k + α · 1|A| ),

where ξ k := [ξk,1 , . . . , ξk,|A| ], and ξk,j denotes the number
of times that action j is played at states belonging to cluster
Ck in the current assignment of a. Explicitly,
X
X X
ξk,j :=
1(at = j) =
1(at = j),
(9)
i∈Ck t:st =i

t:zst =k

which is nothing but the sum of the φi,j ’s of the corresponding states,
X
ξk,j =
φi,j .
i∈Ck

In addition to the actions and control parameters, we now
also need to sample the indicators {zi }, whose conditional
distributions can be expressed in terms of the corresponding
prior model and the likelihood of the triggered actions,
Y
p(zi | z −i , s, a, Θ, α) ∝ p(zi | z −i )
π(at | θ zi ). (10)
t:st =i

2.2.2 Collapsed Gibbs sampling
As before, we derive the collapsed Gibbs sampler by
marginalizing out the control parameters,
Z
p(s, a, z, Θ | α) dΘ
(11)
p(zi | z −i , s, a, α) ∝

∝ p(zi | z −i )
∝ p(zi | z −i )
∝ p(zi | z −i )
∝ p(zi | z −i )

Z

∆K

K
Y

∆K k=1

Z

K
Y

∆K k=1

Z

K
Y

∆K k=1

K
Y

k=1

pθ (θ k | α)
pθ (θ k | α)
pθ (θ k | α)

TY
−1
t=1
|S|

π(at | θ zst ) dΘ

Y Y

i0 =1 t:st =i0

Y

π(at | θ zi0 ) dΘ

Y

i0 :zi0 =k t:st =i0

π(at | θ k ) dΘ



Z
Y
 pθ (θ k | α)
π(at | θ k ) dθ k  .
∆

t:st ∈Ck

Here, we first grouped the actions by their associated states
and then grouped the states themselves by the clusters
{Ck }. Again, this distribution admits an easy sampling
mechanism as it takes the form of a product of Dirichletmultinomials, reweighted by the conditional prior distribution over indicators. In particular, we observe that all
actions played at some state i appear in exactly one of
the K integrals of the last equation. In other words, by
changing the value of zi (i.e. by assigning state i to another
cluster), only two of the involved integrals are affected: the
one belonging to the previously assigned cluster, and the
one of the new cluster. Inference about the value of zi can
thus be carried out using the following two sets of sufficient
statistics:
•

•

φi,j : the number of actions j played at state i,
ψi,j,k : the number of actions j played at states assigned to cluster Ck , excluding state i.

The φi,j ’s are the same as in Eq. (4) and their definition is
repeated here just as a reminder. For the ψi,j,k ’s, on the other
hand, we find the following explicit expression,
X X
ψi,j,k :=
1(at = j),
i0 ∈Ck t:st =i0
i0 6=i

7

which can be also written in terms of the statistics used for
the ordinary Gibbs sampler,

ψi,j,k = ξk,j − 1(i ∈ Ck ) · φi,j .

By collecting these quantities in a vector, i.e. ψ i,k :=
[ψi,1,k , . . . , ψi,|A|,k ], we end up with the following simplified expression,

p(zi = k | z −i , s, a, α) ∝ p(zi = k | z −i ) . . .
K
Y
... ×
D IR M ULT(ψ i,k0 + 1(k 0 = k) · φi | α).
k0 =1

Further, we obtain the following result for the conditional
distribution of action at ,

p(at | a−t , s, z, α) ∝ T (st+1 | st , at ) . . .
Z
Y
π(at0 | θ zst ) dθ zst .
... ×
pθ (θ zst | α)
∆

t0 :zs 0 =zst
t

By introducing the sufficient statistics {ϑt,j }, which count
the number of occurrences of action j among all states that
are currently assigned to the same cluster as st (i.e. the
cluster Czst ), excluding at itself,
X
ϑt,j :=
1(at = j),
t0 :zs 0 =zst
t
t0 6=t

we finally arrive at the following expression,

p(at = j | a−j , s, z, α) ∝ (ϑt,j + α) · T (st+1 | st , at = j).

As for the static model, we can establish a relationship between the statistics used for the ordinary and the collapsed
sampler,
ϑt,j = ξzst ,j − 1(at = j).
2.2.3 Prior models
In order to complete our model, we need to specify a prior
distribution over indicator variables pz (z). The following
paragraphs present three such candidate models:
Non-informative prior
The simplest of all prior models is the non-informative
prior over partitionings, reflecting the assumption that, a
priori, all cluster assignments are equally likely and that the
indicators {zi } are mutually independent. In this case, pz (z)
is constant and, hence, the term p(zi | z −i ) in Eq. (10) and
Eq. (11) disappears, so that the conditional distribution of
indicator zi becomes directly proportional to the likelihood
of the inferred action sequence.
Mixing prior
Another simple yet expressive prior can be realized by the
(finite) Dirichlet mixture model. Instead of assuming that
the indicator variables are independent, the model uses a set
of mixing coefficients q = [q1 , . . . , qK ], where qk represents
the prior probability that an indicator variable takes on
value k . The mixing coefficients are themselves modeled by
a Dirichlet distribution, so that we finally have
γ
· 1K ),
q ∼ D IR(q |
K
zi | q ∼ C AT(zi | q),

where γ is another concentration parameter, controlling the
variability of the mixing coefficients. Note that the indicator
variables are still conditionally independent given the mixing
coefficients in this model. More specifically, for a fixed q , the
conditional distribution of a single indicator in Eq. (10) and
Eq. (11) takes the following simple form,

p(zi = k | z −i , q) = qk .
If the value of q is unknown, we have two options to include
this prior into our model. One is to sample q additionally to
the remaining variables by drawing values from the following conditional distribution during the Gibbs procedure,
|S|
Y
γ
· 1K )
C AT(zi | q)
K
i=1
γ
∝ D IR(q | ζ +
· 1K ),
K
where ζ := [ζ1 , . . . , ζK ], and ζk denotes the number of
variables zi that map to cluster Ck ,

p(q | s, a, z, Θ, α) ∝ D IR(q |

ζk =

|S|
X

1(zi = k).

i=1

Alternatively, we can again make use of the conjugacy property to marginalize out the mixing proportions q during the
inference process, just as we did for the control parameters
in previous sections. The result is (additional) collapsing in
q . In this case, we simply replace the factor p(zi = k | z −i )
in the conditional distribution of zi by
γ
(−i)
(12)
p(zi = k | z −i , γ) ∝ (ζk + ),
K
(−i)

where ζk
is defined like ζk but without counting the
current value of indicator zi ,
(−i)

ζk

:=

|S|
X

i0 =1
i0 6=i

1(zi = k) = ζk − 1(zi = k).

A detailed derivation is omitted here but follows the same
style as for the collapsing in Section 2.1.2.
Spatial prior
Both previous prior models assume (conditional) independence of the indicator variables and, hence, make no specific
assumptions about their dependency structure. However,
we can also use the prior model to promote a certain type of
spatial state clustering. A reasonable choice is, for instance,
to use a model which preferably groups “similar” states
together (in other words, a model which favors clusterings
that assign those states the same local control parameter).
Similarity of states can be expressed, for example, by a
monotonically decreasing decay function f : [0, ∞) → [0, 1]
which takes as input the distance between two states. The
required pairwise distances can be, in turn, defined via some
distance metric χ : S × S → [0, ∞).
In fact, apart from the reasons listed in Section 2.2, there
is an additional motivation, more intrinsically related to the
dynamics of the system, why such a clustering can be useful:
given that the transition model of our system admits locally
smooth dynamics (which is typically the case for real-world

8

systems), the resulting optimal control policy often turns
out to be spatially smooth, too [11]. More specifically, under
an optimal policy, two nearby states are highly likely to
experience similar controls; hence, it is reasonable to assume
a priori that both share the same local control parameter.
For the policy recognition task, it certainly makes sense to
regularize the inference problem by encoding this particular
structure of the solution space into our model. The Potts
model [30], which is a special case of a Markov random
field with pairwise clique potentials [31], offers one way to
do this,

pz (z) ∝

|S|
Y

i=1

exp

!
|S|
βX
f (di,j )δ(zi , zj ) .
2 j=1
j6=i

Here, δ denotes Kronecker’s delta, di,j := χ(si , sj ), i, j ∈
{1, . . . , |S|}, are the state similarity values, and β ∈ [0, ∞)
is the (inverse) temperature of the model which controls
the strength of the prior. From this equation, we can easily
derive the conditional distribution of a single indicator
variable zi as

p(zi | z −i ) ∝ exp β

|S|
X

j=1
j6=i

!

f (di,j )δ(zi , zj ) .

(13)

This completes our inference framework for finite spaces.

2.3

Countably infinite and uncountable state spaces

A major advantage of the clustering approach presented in
the last section is that, due to the limited number of local
policies to be learned from the finite amount of demonstration data, we can now apply the same methodology to
state spaces of arbitrary size, including countably infinite
and uncountable state spaces. This extension had been
practically impossible for the static model because of the
overfitting problem explained in Section 2.2. Nevertheless,
there remains a fundamental conceptual problem: a direct
extension of the model to these spaces would imply that
the distribution over possible state partitionings becomes an
infinite-dimensional object (i.e., in the case of uncountable
state spaces, a distribution over functional mappings from
states to local controllers), requiring an infinite number of
indicator variables. Certainly, such an object is non-trivial to
handle computationally.
However, while the number of latent cluster assignments grows unbounded with the size of the state space,
the amount of observed trajectory data always remains
finite. A possible solution to the problem is, therefore, to
reformulate the inference task on a reduced state space
Se := {s1 , s2 , . . . , sT } containing only states along the observed trajectories. Reducing the state space in this way
means that we need to consider only a finite set of indicator
variables {zt }Tt=1 , one for each expert state s ∈ Se, which
always induces a model of finite size. Assuming that no
state is visited twice, we may further use the same index

set for both variable types.4 In order to limit the complexity
of the dependency structure of the indicator variables for
larger data sets, we further let the value of indicator zt
depend only on a subset of the remaining variables z −t as
defined by some neighborhood rule N . The resulting joint
distribution is then given as
!
T
Y
β X
pz (z | s) ∝
exp
f (dt,t0 )δ(zt , zt0 ) ,
2 t0 ∈N
t=1
t

which now implicitly depends on the state sequence s
through the pairwise distances dt,t0 := χ(st , st0 ), t, t0 ∈
{1, . . . , T } (hence the conditioning on s).
The use of a finite number of indicator variables along
the expert trajectories obviously circumvents the abovementioned problem of representational complexity. Nevertheless, there are some caveats associated with this approach. First of all, using a reduced state space model raises
the question of marginal invariance [32]: if we added a new
trajectory point to the data set, would this change our belief
about the expert policy at previously visited states? In particular, how is this different from modeling that new point
together with the initial ones in the first place? And further,
what does such a reduced model imply for unvisited states?
Can we still use it to make predictions about their local
policies? These questions are, in fact, important if we plan
to use our model to generalize the expert demonstrations to
new situations. For a detailed discussion on this issue, the
reader is referred to the supplement. Here, we focus on the
inferential aspects of the problem, which means to identify
the system parameters at the given trajectory states.
Another (but related) issue resulting from the reduced
modeling approach is that we lose the simple generative
interpretation of the process that could be used to explain
the data generation beforehand. In the case of finite state
spaces, we could think of a trajectory as being constructed
by the following, step-wise mechanism: first, the prior pz (z)
is used to generate a set of indicator variables for all states.
Independently, we pick some value for α from pα (α) and
sample K local control parameters from pθ (θ k | α). To
finally generate a trajectory, we start with an initial state
s1 , generated by p1 (s1 ), select a random action a1 from
π(a1 | s1 , θ zs1 ) and transition to a new state s2 according to
T (s2 | s1 , a1 ), where we select another action a2 , and so on.
Such a directed way of thinking is possible since the finite
model naturally obeys a causal structure where later states
depend on earlier ones and the decisions made there. Furthermore, the cluster assignments and the local controllers
could be generated in advance and isolated from each other
because they were modeled marginally independent.
For the reduced state space model, this interpretation no
longer applies as the model has no natural directionality. In
fact, its variables depend on each other in a cyclic fashion:
altering the value of a particular indicator variable (say, the
one corresponding to the last trajectory point) will have
4. Note that we make this assumption for notational convenience
only and that it is not required from a mathematical point of view.
Nonetheless, for uncountable state spaces the assumption is reasonable
since the event of reaching the same state twice has zero probability
for most dynamic models. In the general case, however, the indicator
variables require their own index set to ensure that each system state is
associated with exactly one cluster, even when visited multiple times.

9

an effect on the values of all remaining indicators due to
their spatial relationship encoded by the “prior distribution”
pz (z | s). Changing the values of the other indicators,
however, will influence the actions being played at the
respective states which, in turn, alters the probability of
ending up with the observed trajectory in the first place and,
hence, the position and value of the indicator variable we
started with. Explaining the data generation of this model
using a simple generative process is, therefore, not possible.
Nevertheless, the individual building blocks of our
model (that is, the policy, the transition model, etc.) together
form a valid distribution over the model variables, which
can be readily used for parameter inference. For the reasons
explained above, it makes sense to define this distribution in
the form of a discriminative model, ignoring the underlying
generative aspects of the process. This is sufficient since we
can always condition on the observed state sequence s,

p(a, Θ, z | s, α) =

K
Y
1
pz (z | s)
pθ (θ k | α) . . .
Zs
k=1

. . . × p1 (s1 )

TY
−1
t=1

T (st+1 | st , at )π(at | θ zst ).

Herein, Zs is a data-dependent normalizing constant. The
structure of this distribution is illustrated by the factor graph
shown in the supplement (Fig. S-1), which highlights the
circular dependence between the variables. Note that, for
any fixed state sequence s, this distribution indeed encodes
the same basic properties as the finite model in Eq. (8).
In particular, the conditional distributions of all remaining
variables remain unchanged, which allows us to apply the
same inference machinery that we already used in the finite
case. For a deeper discussion on the difference between the
two models, we again point to the supplement.

3

N ONPARAMETRIC P OLICY R ECOGNITION

In the last section, we presented a probabilistic policy recognition framework for modeling the expert behavior using
a finite mixture of K local policies. Basically, there are two
situations when such a model is useful:
•
•

either, we know the true number of expert policies,
or, irrespective of the true behavioral complexity, we
want to find an approximate system description in
terms of at most K distinct control situations [27]
(c.f. finite state controllers [29]).

In all other cases, we are faced with the non-trivial problem of choosing K . In fact, the choice of K should not
just be considered a mathematical necessity to perform
inference in our model. By selecting a certain value for
K we can, of course, directly control the complexity class
of potentially inferred expert controllers. However, from a
system identification point of view, it is more reasonable to
infer the required granularity of the state partitioning from
the observed expert behavior itself, instead of enforcing a
particular model complexity. This way, we can gain valuable
information about the underlying control structure and state
representation used by the expert, which offers a possibility
to learn a state partitioning of task-appropriate complexity
directly from the demonstration data. Hence, the problem

of selecting the right model structure should be considered
as part of the inference problem itself.
From a statistical modeling perspective, there are two
common ways to approach this problem. One is to make use
of model selection techniques in order to determine the most
parsimonious model that is in agreement with the observed
data. However, choosing a particular model complexity still
means that we consider only one possible explanation for
the data, although other explanations might be likewise
plausible. For many inference tasks, including this one, the
more elegant approach is to keep the complexity flexible
and, hence, adaptable to the data. Mathematically, this can
be achieved by assuming a potentially infinite set of model
parameters (in our case controllers) from which we activate
only a finite subset to explain the particular data set at hand.
This alternative way of thinking opens the door to the rich
class of nonparametric models, which provide an integrated
framework to formulate the inference problem over both
model parameters and model complexity as a joint learning
problem.
3.1

A Dirichlet process mixture model

The classical way to nonparametric clustering is to use
a Dirichlet process mixture model (DPMM) [33]. These
models can be obtained by starting from a finite mixture
model and letting the number of mixture components (i.e.
the number of local controllers) approach infinity. In our
case, we start with the clustering model from Section 2.2,
using a mixing prior over indicator variables,

q ∼ D IR(q |

γ
K

· 1K )

θ k ∼ D IR(θ k | α · 1|A| )
s1 ∼ p1 (s1 )

zi | q ∼ C AT(zi | q)

at | st , Θ, z ∼ π(at | θ zst ) (14)
st+1 | st , at ∼ T (st+1 | st , at ) .

From these equations, we arrive at the corresponding nonparametric model as K goes to infinity. For the theoretical
foundations of this limit, the reader is referred to the more
general literature on Dirichlet processes, such as [33], [34].
In this paper, we restrict ourselves to providing the resulting
sampling mechanisms for the policy recognition problem.
In a DPMM, the mixing proportions q of the local
parameters are marginalized out (that is, we use a collapsed sampler). The resulting distribution over partitionings is described by a Chinese restaurant process (CRP) [35]
which can be derived, for instance, by considering the limit
K → ∞ of the mixing process induced by the Gibbs update
in Eq. (12),
(
(−i)
ζ
if k ∈ {1, . . . , K ∗ },
p(zi = k | z −i , γ) ∝ k
(15)
γ
if k = K ∗ + 1.
Here, K ∗ denotes the number of distinct entries in z −i
which are represented by the numerical values {1, . . . , K ∗ }.
In this model, a state joins an existing cluster (i.e. a group of
states whose indicators have the same value) with probability proportional to the number of states already contained in
that cluster. Alternatively, it may create a new cluster with
probability proportional to γ .
From the model equations (14) it is evident that, given a
particular setting of indicators, the conditional distributions
of all other variable types remain unchanged. Effectively, we

10

only replaced the prior model pz (z) by the CRP. Hence, we
can apply the same Gibbs updates for the actions and controllers as before and need to rederive only the conditional
distributions of the indicator variables under consideration
of the above defined process. According to Eq. (15), we
herein need to distinguish whether an indicator variable
takes a value already occupied by other indicators (i.e. it
joins an existing cluster) or it is assigned a new value (i.e.
∗
it creates a new cluster). Let {θ k }K
k=1 denote the set of
control parameters associated with z −i . In the first case
(k ∈ {1, . . . , K ∗ }), we can then write

S

c2

s2

c8

s3

s4

c3
s1

c1

s7
c4

c6
s6

s5
c5

s8
c7

c9
s9

Fig. 3: Schematic illustration of the ddCRP-based clustering
applied to the reduced state space model in Section 2.3.
p(zi = k |
Each trajectory state is connected to some other state of
= p(zi = k | z −i , {at }t:st =i , θ k , α, γ)
the sequence. The connected components of the resulting
∝ p(zi = k | z −i , θ k , α, γ)p({at }t:st =i | zi = k, z −i , θ k , α, γ) graph implicitly define the state clustering. Coloring of the
∝ p(zi = k | z −i , γ)p({at }t:st =i | θ k )
background illustrates the spatial cluster extrapolation (see
Y
(−i)
Section A in the supplement). Note that the underlying
∝ ζk ·
π(at | θ k ).
decision-making process is assumed to be discrete in time;
t:st =i
the continuous gray line shown in the figure is only to
In the second case (k = K ∗ + 1), we instead obtain
indicate the temporal ordering of the trajectory states.
∗
K∗
p(zi = K + 1 | z −i , s, a, {θ k }k=1 , α, γ)
∗
z −i , s, a, {θ k0 }K
k0 =1 , α, γ)

= p(zi = K ∗ + 1 | z −i , {at }t:st =i , α, γ)
∝ p(zi = K ∗ + 1 | z −i , α, γ) . . .
. . . × p({at }t:st =i | zi = K ∗ + 1, z −i , α, γ)
∝ p(zi = K ∗ + 1 | z −i , γ)p({at }t:st =i | zi = K ∗ + 1, α)
Z
∝γ·
p({at }t:st =i | θ K ∗ +1 )pθ (θ K ∗ +1 | α) dθ K ∗ +1
Z∆ Y
π(at | θ K ∗ +1 )pθ (θ K ∗ +1 | α) dθ K ∗ +1
∝γ·
∆ t:s =i
t

∝ γ · D IR M ULT(φi | α).

If a new cluster is created, we further need to initialize the
corresponding control parameter θ K ∗ +1 by performing the
respective Gibbs update, i.e. by sampling from
∗

p(θ K ∗ +1 | z, s, a, {θ k }K
k=1 , α, γ)
= p(θ K ∗ +1 | {at }t:zst =K ∗ +1 , α)
∝ pθ (θ K ∗ +1 | α)p({at }t:zst =K ∗ +1 | θ K ∗ +1 )
Y
∝ pθ (θ K ∗ +1 | α) ·
π(at | θ K ∗ +1 )
t:zst =K ∗ +1

∝ D IR(θ K ∗ +1 | ξ K ∗ +1 + α · 1|A| ).

Should a cluster get unoccupied during the sampling process, the corresponding control parameter may be removed
from the stored parameter set {θ k } and the index set for
k needs to be updated accordingly. Note that this sampling
mechanism is a specific instance of Algorithm 2 described in
[33]. A collapsed variant can be derived in a similar fashion.
3.2 Policy recognition using the distance-dependent
Chinese restaurant process
In the previous section, we have seen that the DPMM can
be derived as the nonparametric limit model of a finite
mixture using a set of latent mixing proportions q for the
clusters. Although the DPMM allows us to keep the number
of active controllers flexible and, hence, adaptable to the
complexity of the demonstration data, the CRP as the underlying clustering mechanism does not capture any spatial

dependencies between the indicator variables. In fact, in the
CRP, the indicators {zi } are coupled only via their relative
frequencies (see Eq. (15)) but not through their individual
locations in space, resulting in an exchangeable collection
of random variables [35]. In fact, one could argue that the
spatial structure of the clustering problem is a priori ignored.
The fact that DPMMs are nevertheless used for spatial
clustering tasks can be explained by the particular form
of data likelihood models that are typically used for the
mixture components. In a Gaussian mixture model [36], for
instance, the spatial clusters emerge due to the unimodal
nature of the mixture components, which encodes the locality property of the model needed to obtain a meaningful
spatial clustering of the data. For the policy recognition
problem, however, the DPMM is not able to exploit any
spatial information via the data likelihood since the clustering of states is performed at the level of the inferred action
information (see Eq. (10)) and not on the state sequence
itself. Consequently, we cannot expect to obtain a smooth
clustering of the system state space, especially when the
expert policies are overlapping (i.e. when they share one or
more common actions) so that the action information alone
is not sufficient to discriminate between policies. For uncountable state spaces, this problem is further complicated
by the fact that we observe at most one expert state transition
per system state. Here, the spatial context of the data is the
only information which can resolve this ambiguity.
In order to facilitate a spatially smooth clustering, we
therefore need to consider non-exchangeable distributions
over partitionings. More specifically, we need to design our
model in such a way that, whenever a state s is “close” to
some other state s0 and assigned to some cluster Ck , then,
a priori, s0 should belong to the same cluster Ck with high
probability. In that sense, we are looking for the nonparametric counterpart of the Potts model. One model with such
properties is the distance-dependent Chinese restaurant

11

process (ddCRP) [32].5 As opposed to the traditional CRP,
the ddCRP explicitly takes into account the spatial structure
of the data. This is done in the form of pairwise distances
between states, which can be obtained, for instance, by
defining an appropriate distance metric on the state space
(see Section 2.2.3). Instead of assigning states to clusters as
done by the CRP, the ddCRP assigns states to other states
according to their pairwise distances. More specifically, the
probability that state i gets assigned to state j is defined as
(
ν
if i = j,
p(ci = j | D, ν) ∝
(16)
f (di,j ) otherwise,
where ν ∈ [0, ∞) is called the self-link parameter of the process, D denotes the collection of all pairwise state distances,
and ci is the “to-state” assignment of state i, which can be
thought of as a directed edge on the graph defined on the
set of all states (see Fig. 3). Accordingly, i and j in Eq. (16)
can take values {1, . . . , |S|} for the finite state space model
and {1, . . . , T } for our reduced state space model. The state
clustering is then obtained as a byproduct of this mapping
via the connected components of the resulting graph (see
Fig. 3 again).
Replacing the CRP by the ddCRP and following the
same line of argument as in [32], we obtain the required
conditional distribution of the state assignment ci as


j = i,
ν
p(ci = j | c−i ,s,a, α,D, ν) ∝ f (di,j )
no clusters merged,


f (di,j )·L Czi and Czj merged,

where we use the shorthand notation

L=

D IR M ULT(ξ zi + ξ zj | α)

D IR M ULT(ξ zi | α)D IR M ULT(ξ zj | α)

for the data likelihood term. The ξk,j ’s are defined as in
Eq. (9) but are based on the clustering which arises when we
ignore the current link ci . The resulting Gibbs sampler is a
collapsed one as the local control parameters are necessarily
marginalized out during the inference process.

4

S IMULATION R ESULTS

In this section, we present simulation results for two types of
system dynamics. As a proof-of-concept, we first investigate
the case of uncountable state spaces which we consider the
more challenging setting for reasons explained earlier. To
compare our framework with existing methods, we further
provide simulation results for the standard grid world
benchmark (see e.g. [9], [11], [19]). It should be pointed
out, however, that establishing a fair comparison between
LfD models is generally difficult due to their different
working principles (e.g. reward prediction vs. action prediction), objectives (system identification vs. optimal control),
requirements (e.g. MDP solver, knowledge of the expert’s
discount factor, countable vs. uncountable state space), and
5. Note that the authors of [32] avoid calling this model nonparametric since it cannot be cast as a mixture model originating from a random
measure. However, we stick to this term in order to make a clear
distinction to the parametric models in Section 2, and to highlight the
fact that there is no parameter K determining the number of controllers.

Fig. 4: Schematic illustration of the expert policy used in
Section 4.1, which applies eight local controllers to sixteen
distinct regions. A sample trajectory is shown in color.
assumptions (e.g. deterministic vs. stochastic expert behavior). Accordingly, our goal is rather to demonstrate the
prediction abilities of the considered models than to push
the models to their individual limits. Therefore, and to reduce the overall computational load, we tuned most model
hyper-parameters by hand. Our code is available at https:
//github.com/AdrianSosic/BayesianPolicyRecognition.
4.1

Example 1: uncountable state space

As an illustrative example, we consider a dynamical system which describes the circular motion of an agent on
a two-dimensional state space. The actions of the agent
correspond to 24 directions that divide the space of possible
angles [0, 2π) into equally-sized intervals. More specifically,
action j corresponds to the angle (j − 1) 2π
24 . The transition
model of the system is defined as follows: for each selected
action, the agent first makes a step of length µ = 1 in
the intended direction. The so-obtained position is then
distorted by additive zero-mean isotropic Gaussian noise of
variance σ 2 . This defines our transition kernel as

T (st+1 | st , at = j) = N (st+1 | st + µ · ej , σ 2 I),

(17)

where st , st+1 ∈ R2 , ej denotes the two-dimensional unit
vector pointing in the direction of action j , and I is the
two-dimensional identity matrix. The overall goal of our
agent is to describe a circular motion around the origin in
the best possible manner allowed by the available actions.
However, due to limited sensory information, the agent is
not able to observe its exact position on the plane but can
only distinguish between certain regions of the state space,
as illustrated by Fig. 4. Also, the agent is unsure about the
optimal control strategy, i.e. it does not always make optimal
decisions but selects its actions uniformly at random from a
subset of actions, consisting of the optimal one and the two
actions pointing to neighboring directions (see Fig. 4 again).
To increase the difficulty of the prediction task, we further
let the agent change the direction of travel whenever the
critical distance of r = 5 to the origin is exceeded.
Having defined the expert behavior, we generate 10
sample trajectories of length T = 100. Herein, we assume
a motion noise level of σ = 0.2 and initialize the agent’s

12

0

0

χ(s, s ) = ||s − s ||2 .

(18)

The corresponding similarity values are calculated using a
Gaussian-shaped kernel. More specifically,
!
d2
fPotts (d) = exp − 2
σf
for the Potts model and

fddCRP (d) = (1 − )fPotts (d) + 

for the ddCRP model, with σf = 1 and a constant offset of
 = 0.01 which ensures that states with large distances can
still join the same cluster. For the Potts model, we further
use a neighborhood structure containing the eight closest
trajectory points of a state. This way, we ensure that, in
principle, each local expert policy may occur at least once
in the neighborhood of a state. The concentration parameter
for the local controls is set to α = 1, corresponding to a
uniform prior belief over local policies.
A major drawback of the Potts model is that posterior
inference about the temperature parameter β is complicated
due to the nonlinear effect of the parameter on the normalization of the model. Therefore, we manually selected
a temperature of β = 1.6 based on a minimization of
the average policy prediction error (discussed below) via
parameter sweeping. As opposed to this, we extend the
inference problem for the ddCRP to the self-link parameter ν
as suggested in [32]. For this, we use an exponential prior,

pν (ν) = E XP(ν | λ),

with rate parameter λ = 0.1, and applied the independence
Metropolis-Hastings algorithm [37] using pν (ν) as proposal
distribution with an initial value of ν = 1. In all our
simulations, the sampler quickly converged to its stationary
distribution, yielding posterior values for ν with a mean of
0.024 and a standard deviation of 0.023.
To locally compare the predicted policy with the ground
truth at a given state, we compute their earth mover’s
distance (EMD) [38] with a ground distance metric measuring the absolute angular difference between the involved
actions. To track the learning progress of the algorithms,
we calculate the average EMD over all states of the given
trajectory set at each Gibbs iteration. Herein, the local policy
predictions are computed from the single Gibbs sample of
the respective iteration, consisting of all sampled actions, indicators and – in case of non-collapsed sampling – the local
control parameters. The resulting mean EMDs and standard
deviations are depicted in Fig. 5. The inset further shows
the average EMD computed at non-trajectory states which
are sampled on a regular grid (depicted in the supplement),
reflecting the quality of the resulting spatial prediction.
As expected, the finite mixture model (using the true
number of local policies, a collapsed mixing prior, and

1.4
1.4
1.2
1
0.8
0.6
0.4
10 0

1.2

average EMD

position uniformly at random on the unit circle. An example
trajectory is shown in Fig. 4. The obtained trajectory data is
fed into the presented inference algorithms to approximate
the posterior distribution over expert controllers, and the
whole experiment is repeated in 100 Monte Carlo runs.
For the spatial models, we use the Euclidean metric to
compute the pairwise distances between states,

1
0.8
0.6
0.4
10 0

10 1

10 2

mixing prior
Potts (complete)
Potts (collapsed)
ddCRP

10 1

10 2

Gibbs iteration

Fig. 5: Average policy prediction error at the simulated
trajectory states (main figure) and at non-trajectory states
(inset). Shown are the empirical mean values and standard
deviations, estimated from 100 Monte Carlo runs.

γ = 1) is not able to learn a reasonable policy representation
from the expert demonstrations since it does not explore the
spatial structure of the data. In fact, the resulting prediction
error shows only a slight improvement as compared to
an untrained model. In contrast to this, all spatial models
capture the expert behavior reasonably well. In agreement
with our reasoning in Section 2.1.2, we observe that the
collapsed Potts model mixes significantly faster and has a
smaller prediction variance than the non-collapsed version.
However, the ddCRP model gives the best result, both
in terms of mixing speed (see [32] for an explanation of
this phenomenon) and model accuracy. Interestingly, this is
despite the fact that the ddCRP model additionally needs to
infer the number of local controllers necessary to reproduce
the expert behavior. The corresponding posterior distribution, which shows a pronounced peak at the true number,
is depicted in the supplement. There, we also provide additional simulation results which give insights into the learned
state partitioning and the resulting spatial policy prediction
error. The results reveal that all expert motion patterns can
be identified by our algorithm.
4.2

Example 2: finite state space

In this section, we compare the prediction capabilities of our
model to existing LfD frameworks, in particular: the maximum margin method in [9] (max-margin), the maximum
entropy approach in [12] (max-ent), and the expectationmaximization algorithm in [11] (EM). For the comparison,
we restrict ourselves to the ddCRP model which showed
the best performance among all presented models.
As a first experiment, we compare all methods on a finite
version of the setting in Section 4.1, which is obtained by
discretizing the continuous state space into a regular grid
S = {(x, y) ∈ Z2 : |x|, |y| ≤ 10}, resulting in a total of 441
states. The transition probabilities are chosen proportional
to the normal densities in Eq. (17) sampled at the grid points.
Here, we used a noise level of σ = 1 and a reduced number
of eight actions. Probability mass “lying outside” the finite
grid area is shifted to the closest border states of the grid.
Figure 6a delineates the average EMD over the number
of trajectories (each of length T = 10) provided for training.

13

1.5

1.5

average EMD

average EMD

2

1
0.5
0
10 1

max-ent
max-margin
EM
ddCRP

10 2

10 3

10 4

number of trajectories

(a) learning curves (circular policy)

1.5

max-ent
max-margin
EM
ddCRP

1

1

0.5

0
10 1

max-ent
max-margin
EM
ddCRP

0.5

10 2

10 3

10 4

0
0

0.2

0.4

0.6

0.8

1

number of trajectories

(b) learning curves (MDP policy)

(c) model robustness (MDP policy)

Fig. 6: Average EMD values for the prediction task described in Section 4.2: (a) circular policy (b,c) MDP policy. Shown are
the empirical mean values and standard deviations, estimated from 100 Monte Carlo runs. The EMD values are computed
based on (a) the predicted action distributions and (b,c) the predicted next-state distributions. Note that the curves of
max-ent (purple) and max-margin (yellow) in subfigure (a) lie on top of each other.
We observe that neither of the two intentional models (maxent and max-margin) is able to capture the demonstrated expert behavior. This is due to the fact that the circular expert
motion cannot be explained by a simple state-dependent
reward structure but requires a more complex state-action
reward model, which is not considered in the original formulations [9], [12]. While the EM model is indeed able to
capture the general trend of the data, the prediction is less
accurate as compared to that of the ddCRP model, since it
cannot reproduce the stochastic nature of the expert policy.
In fact, this difference in performance will become even
more pronounced for expert policies which distribute their
probability mass on a larger subset of actions. Therefore,
the ddCRP model outperforms all other models since the
provided expert behavior violates their assumptions.
To analyze how the ddCRP competes against the other
models in their nominal situations, we further compare all
algorithms on a standard grid world task where the expert
behavior is obtained as the optimal response to a simple
state-dependent reward function. Herein, each state on the
grid is assigned a reward with a chance of 1%, which is then
drawn from a standard normal distribution. Worlds which
contain no reward are discarded. The discount factor 0.9,
which is used to compute the expert policy (see [17]), is
provided as additional input for the intentional models.
The results are shown in Figure 6b, which illustrates that
the intention-based max-margin method outperforms all
other methods for small amounts of training data. The subintentional methods (EM and ddCRP), on the other hand,
yield better asymptotic estimates and smaller prediction
variances. It should be pointed out that the three reference
methods have a clear advantage over the ddCRP in this case
because they assume a deterministic expert behavior a priori
and do not need to infer this piece of information from the
data. Despite this additional challenge, the ddCRP model
yields a competitive performance.
Finally, we compare all approaches in terms of their
robustness against modeling errors. For this purpose, we
repeat the previous experiment with a fixed number of 1000
trajectories but employ a different transition model for inference than used for data generation. More specifically, we utilize an overly fine-grained model consisting of 24 directions,
assuming that the true action set is unknown, as suggested

in Section 1.1. Additionally, we perturb the assumed model
by multiplying (and later renormalizing) each transition
probability with a random number generated according to
f (u) = tan( π4 (u + 1)), with u ∼ U NIFORM(−η, η) and
perturbation strength η ∈ [0, 1]. Due the resulting model
mismatch, a comparison to the ground truth policy based
on the predicted action distribution becomes meaningless.
Instead, we compute the Euclidean EMDs between the
true and the predicted next-state distributions, which we
obtain by marginalizing the actions of the true/assumed
transition model with respect to the true/learned policy.
Figure 6c depicts the resulting prediction performance for
different perturbation strengths η . The results confirm that
our approach is not only less sensitive to modeling errors as
argued in Section 1.1; also, the prediction variance is notably
smaller than those of the intentional models.

5

C ONCLUSION

In this work, we proposed a novel approach to LfD by
jointly learning the latent control policy of an observed
expert demonstrator together with a task-appropriate representation of the system state space. With the described
parametric and nonparametric models, we presented two
formulations of the same problem that can be used either
to learn a global system controller of specified complexity,
or to infer the required model complexity from the observed expert behavior itself. Simulation results for both
countable and uncountable state spaces and a comparison
to existing frameworks demonstrated the efficacy of our
approach. Most notably, the results showed that our method
is able to learn accurate predictive behavioral models in
situations where intentional methods fail, i.e. when the
expert behavior cannot be explained as the result of a simple
planning procedure. This makes our method applicable to
a broader range of problems and suggests its use in a more
general system identification context where we have no such
prior knowledge about the expert behavior. Additionally,
the task-adapted state representation learned through our
framework can be used for further reasoning.

14

R EFERENCES
[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]

B. D. Argall, S. Chernova, M. Veloso, and B. Browning, “A survey
of robot learning from demonstration,” Robotics and Autonomous
Systems, vol. 57, no. 5, pp. 469–483, 2009.
M. P. Deisenroth, D. Fox, and C. E. Rasmussen, “Gaussian processes for data-efficient learning in robotics and control,” IEEE
Transactions on Pattern Analysis and Machine Intelligence, vol. 37,
no. 2, pp. 408–423, 2015.
P. Abbeel and A. Y. Ng, “Exploration and apprenticeship learning
in reinforcement learning,” in Proc. 22nd International Conference on
Machine Learning, 2005, pp. 1–8.
D. Michie, M. Bain, and J. Hayes-Miches, “Cognitive models from
subcognitive skills,” IEE Control Engineering Series, vol. 44, pp. 71–
99, 1990.
C. Sammut, S. Hurst, D. Kedzier, and D. Michie, “Learning to fly,”
in Proc. 9th International Workshop on Machine Learning, 1992, pp.
385–393.
P. Abbeel, A. Coates, and A. Y. Ng, “Autonomous helicopter aerobatics through apprenticeship learning,” The International Journal
of Robotics Research, 2010.
D. A. Pomerleau, “Efficient training of artificial neural networks
for autonomous navigation,” Neural Computation, vol. 3, no. 1, pp.
88–97, 1991.
C. G. Atkeson and S. Schaal, “Robot learning from demonstration,” in Proc. 14th International Conference on Machine Learning,
vol. 97, 1997, pp. 12–20.
P. Abbeel and A. Y. Ng, “Apprenticeship learning via inverse
reinforcement learning,” in Proc. 21st International Conference on
Machine Learning, 2004.
A. Panella and P. J. Gmytrasiewicz, “Nonparametric Bayesian
learning of other agents’ policies in interactive POMDPs,” in
Proc. International Conference on Autonomous Agents and Multiagent
Systems, 2015, pp. 1875–1876.
A. Šošić, A. M. Zoubir, and H. Koeppl, “Policy recognition via expectation maximization,” in Proc. 41st IEEE International Conference
on Acoustics, Speech and Signal Processing, 2016.
B. D. Ziebart, A. L. Maas, J. A. Bagnell, and A. K. Dey, “Maximum entropy inverse reinforcement learning.” in Proc. 23rd AAAI
Conference on Artificial Intelligence, 2008, pp. 1433–1438.
A. Y. Ng and S. J. Russell, “Algorithms for inverse reinforcement
learning.” in Proc. 17th International Conference on Machine Learning,
2000, pp. 663–670.
D. Ramachandran and E. Amir, “Bayesian inverse reinforcement
learning,” Proc. 20th International Joint Conference on Artical Intelligence, vol. 51, pp. 2586–2591, 2007.
S. D. Parsons, P. Gymtrasiewicz, and M. J. Wooldridge, Game theory
and decision theory in agent-based systems.
Springer Science &
Business Media, 2012, vol. 5.
K. Hindriks and D. Tykhonov, “Opponent modelling in automated
multi-issue negotiation using Bayesian learning,” in Proc. 7th
International Joint Conference on Autonomous Agents and Multiagent
Systems, 2008, pp. 331–338.
R. S. Sutton and A. G. Barto, Reinforcement learning: An introduction.
MIT press, 1998.
U. Syed and R. E. Schapire, “A game-theoretic approach to apprenticeship learning,” in Advances in Neural Information Processing
Systems, 2007, pp. 1449–1456.
B. Michini and J. P. How, “Bayesian nonparametric inverse reinforcement learning,” in Machine Learning and Knowledge Discovery
in Databases. Springer, 2012, pp. 148–163.
C. G. Atkeson and J. C. Santamaria, “A comparison of direct and
model-based reinforcement learning,” in International Conference
on Robotics and Automation, 1997.
C. A. Rothkopf and C. Dimitrakakis, “Preference elicitation and inverse reinforcement learning,” in Machine Learning and Knowledge
Discovery in Databases. Springer, 2011, pp. 34–48.
O. Pietquin, “Inverse reinforcement learning for interactive systems,” in Proc. 2nd Workshop on Machine Learning for Interactive
Systems, 2013, pp. 71–75.
S. Schaal, “Is imitation learning the route to humanoid robots?”
Trends in Cognitive Sciences, vol. 3, no. 6, pp. 233–242, 1999.
K. Dvijotham and E. Todorov, “Inverse optimal control with
linearly-solvable MDPs,” in Proc. 27th International Conference on
Machine Learning, 2010, pp. 335–342.
E. Charniak and R. P. Goldman, “A Bayesian model of plan
recognition,” Artificial Intelligence, vol. 64, no. 1, pp. 53–79, 1993.

[26] B. Piot, M. Geist, and O. Pietquin, “Learning from demonstrations:
Is it worth estimating a reward function?” in Machine Learning and
Knowledge Discovery in Databases. Springer, 2013, pp. 17–32.
[27] M. Waltz and K. Fu, “A heuristic approach to reinforcement
learning control systems,” IEEE Transactions on Automatic Control,
vol. 10, no. 4, pp. 390–398, 1965.
[28] F. Doshi-Velez, D. Pfau, F. Wood, and N. Roy, “Bayesian nonparametric methods for partially-observable reinforcement learning,” IEEE Transactions on Pattern Analysis and Machine Intelligence,
vol. 37, no. 2, pp. 394–407, 2015.
[29] N. Meuleau, L. Peshkin, K.-E. Kim, and L. P. Kaelbling, “Learning
finite-state controllers for partially observable environments,” in
Proc. 15th Conference on Uncertainty in Artificial Intelligence, 1999,
pp. 427–436.
[30] R. B. Potts, “Some generalized order-disorder transformations,”
in Mathematical Proceedings of the Cambridge Philosophical Society,
vol. 48, no. 01, 1952, pp. 106–109.
[31] D. Koller and N. Friedman, Probabilistic Graphical Models: Principles
and Techniques. MIT press, 2009.
[32] D. M. Blei and P. I. Frazier, “Distance dependent Chinese restaurant processes,” The Journal of Machine Learning Research, vol. 12,
pp. 2461–2488, 2011.
[33] R. M. Neal, “Markov chain sampling methods for Dirichlet process
mixture models,” Journal of Computational and Graphical Statistics,
vol. 9, no. 2, pp. 249–265, 2000.
[34] T. S. Ferguson, “A Bayesian analysis of some nonparametric
problems,” The Annals of Statistics, pp. 209–230, 1973.
[35] D. J. Aldous, Exchangeability and related topics. Springer, 1985.
[36] C. E. Rasmussen, “The infinite Gaussian mixture model,” in Advances in Neural Information Processing Systems. MIT Press, 2000,
pp. 554–560.
[37] S. Chib and E. Greenberg, “Understanding the MetropolisHastings algorithm,” The American Statistician, vol. 49, no. 4, pp.
327–335, 1995.
[38] Y. Rubner, C. Tomasi, and L. J. Guibas, “A metric for distributions
with applications to image databases,” in 6th International Conference on Computer Vision, 1998, 1998, pp. 59–66.

Adrian Šošić is a member of the Signal Processing Group and an associate member of
the Bioinspired Communication Systems Lab at
Technische Universität Darmstadt. Currently, he
is working towards his Ph.D. degree under the
supervision of Prof. Abdelhak M. Zoubir and
Prof. Heinz Koeppl. His research interests center
around topics from machine learning and (inverse) reinforcement learning, with a focus on
probabilistic inference, multi-agent systems, and
Bayesian nonparametrics.
Abdelhak M. Zoubir is professor at the Department of Electrical Engineering and Information Technology at Technische Universität Darmstadt, Germany. His research interest lies in statistical methods for signal processing with emphasis on bootstrap techniques, robust detection
and estimation, and array processing applied
to telecommunications, radar, sonar, automotive
monitoring and biomedicine.

Heinz Koeppl is professor at the Department of
Electrical Engineering and Information Technology at Technische Universität Darmstadt, Germany. His research interests include Bayesian
inference methods for biomolecular data and
methods for reconstructing large-scale biological
or technological multi-agent systems from observational data.

1

A Bayesian Approach to Policy Recognition
and State Representation Learning
(Supplement)
Adrian Šošić, Abdelhak M. Zoubir and Heinz Koeppl

F

A
IN

M ARGINAL I NVARIANCE & P OLICY P REDICTION
L ARGE S TATE S PACES

z2

When we extended our reasoning to large state spaces
in Section 2.3 using a reduced state space model (see
Fig. S-1), we inevitably arrived at the following questions:
By modeling the expert behavior only along observed
trajectories, what does the resulting model imply for the
remaining states of the state space? Can we still use it for
predicting their local policies? The purpose of this section
is to provide an in-depth discussion on the implications of
this reduced modeling approach in the context of policy
prediction.
When investigating the above-mentioned questions from a
probabilistic perspective (i.e. by analyzing the induced joint
distribution of our model), it turns out that they are strongly
related to what is known as marginal invariance [1] (sometimes also referred to as marginalization property or simply
consistency [2]). This property states that a model is consistent in the sense that it always provides the same marginal
distributions for any subset of its variables, irrespective of
the initial model size. In other words, a marginally invariant
policy model yields the same answer for the given trajectory
points, even if we include additional states into our reduced
set Se for which we have not observed any demonstrations.
For our spatial models, that is, the Potts model and
the ddCRP, it can be shown that this consistency property
is indeed lacking (see [1] for a detailed discussion). This
means that we cannot expect to get compatible results when
conducting our reduced model inference on two data sets
of different sizes. On the contrary, making predictions for
new states would always require to rerun our Gibbs sampler
on the augmented data set, including all additional states.
This brings us to the following practical dilemma: imagine
an on-line policy recognition scenario where we observe an
expert controlling our system. After a certain period of time,
we are asked to take over control, using the experience we
have acquired during the observation period. Each control
command, whether performed by the expert or by us, will
trigger a new state transition, meaning that new data points
arrive sequentially one after another. Consequently, it is
impossible to decide in advance which states to include in

S

s2
s3
z3

s1

z4

s7

s4

z1
z5

s5

s8
z8

z7
s6
z6

z9
s9

(a) reduced state space model

zt

zt+1

zt+2

pz

st

st+1

T

at

at+1

st+2

T

at+2

π
π

π
α
pθ

Θ

(b) corresponding factor graph

Fig. S-1: (a) Illustration of the reduced state space model,
which operates on the space Se = {s1 , s2 , . . . , sT } of visited
trajectory states. Note that the underlying decision-making
process is assumed to be discrete in time; the continuous
gray line shown in the figure is only to indicate the temporal
ordering of the trajectory states. (b) Corresponding factor
graph, highlighting the circular dependence between the
variables. The factors are defined by the same building
blocks that are used for the finite state space model. Observed variables are shaded in gray.

2
-15

1.6

-10

1.2
0
1
5

5

-5

0

5

10

-10

7

15

-15 -10

(a) mean prediction error

6
5

0
4
5

3

10

2

15

1

0.2

15

0.4
-15 -10

0.4

10

0.6

15

0.6

0

0.8

10

8

-5

-5
EMD

-5

-15

0.8

1.4

EMD

-10

1

controller

1.8

-15

-5

0

5

10

15

-15 -10

(b) standard deviation

-5

0

5

10

15

(c) sample partitioning

0.4

probability

probability

0.6

0.3
0.2
0.1

5

10

15

0.2
0
0

0
0

0.4

20

5

number of controllers

10
action 15
s

(d) inferred model complexity

20
25

8
7
6 rs
5
4 rolle
3
2 cont
1

(e) posterior sample of the local controllers

Fig. S-2: Simulation results for the ddCRP model on the continuous state space task described in Section 4.1. (a) Mean
values of the spatial policy prediction error. (b) Standard deviations of the spatial policy prediction error. (c) Example
partitioning of the state space, based on the local controllers depicted in sub-figure (e). The expert partitioning is shown
in Fig. 4 of the main paper. (d) Posterior distribution of the number of local controllers. (e) Posterior sample of the local
controllers found by the model. The results in (a,b,d) are based on 100 Monte Carlo runs while (c,e) are obtained from a
single posterior sample. The figures in the top row were rendered using a spatial resolution of 2000x2000.
our reduced space Se and which not. A rigorous approach in
the above-described sense would thus require to recalibrate
the model after each state transition – a costly operation.
However, it is evident that the resulting data set is
naturally divided into two disjoint parts, namely the expert
demonstrations and the subsequent states reached during
execution of the learned policy. Clearly, transitions occurring
after the learning phase should by no means affect our
belief about the expert policy and, hence, they should be
completely discarded from the model. The easiest way to
achieve this is, indeed, to “freeze” the model after the
demonstration phase and to use the learned parameters to
extrapolate the gathered policy information to surrounding
states. This can be done, for instance, by retaining the
structure of the involved spatial prior model to compute the
resulting maximum a posteriori estimates for the extrapolated
indicators of the new states, based on the inferred model
parameters. In the case of the ddCRP, this coincides with
the nearest-neighbor estimate (see Eq. (16)),

ĉnew = arg max f (dt,new ) = arg min dt,new .
t∈{1,...,T }

(S-1)

t∈{1,...,T }

Herein, ĉnew is the estimate for the indicator of the new
state and dt,new denotes the distance of that state to the tth
trajectory point.

Now, one could argue that the comfort of retaining a
finite model structure for modeling inference problems on
countably infinite or uncountable state spaces comes at the
cost of not being able to provide a consistent posterior
predictive distribution. However, the reduced state space
approach allows us to incorporate the spatial information
of the data in a fairly natural manner (i.e. in the form of
pairwise distances), providing an easy way to model the
expert behavior. Furthermore, our results demonstrate that
the reduced model is able to capture the relevant spatial
properties of a policy sufficiently accurate in order to make
profound predictions about unseen states (see also subsequent section). Whether there exist alternative tractable
models with similar properties remains to be seen.

B

A DDITIONAL S IMULATION R ESULTS

In this section, we provide additional simulation results
for the ddCRP model on the continuous state space task
described in Section 4.1.
Figure S-2a visualizes the spatial EMD prediction errors
of the trained model in the form of a heat map, which
compares the ground truth expert policy at non-trajectory
points with the mean prediction provided by our model.
The test points are placed on a regular grid of size 2000x2000
centered around the origin. The required indicator variables

3

at the interpolated states are computed according to Eq. (S1). In line with our expectation, the prediction error reaches
its maximum at the policy boundaries but is comparably
small within each policy region, indicating a good model
fit. Note that the “windmill shape” of the error can be
explained as a result of the reduced state space approach in
combination with the inherent asymmetry of the used data
generation scheme: regions of the state space containing trajectory endings are locally underrepresented in the data set
(see example trajectory in Fig. 4 in the paper); this increases
the chance of assigning the end points of a trajectory to the
cluster of the preceding region, resulting in a smearing of
the previous cluster into the next region.
Also, we can observe that the variance of the error
(Fig. S-2b) reaches its maximum at the transition regions
and generally grows with the distance to the supporting
trajectory data, reflecting the increasing prediction uncertainty at cluster boundaries and regions far from the expert
demonstrations. Both figures were computed based on the
learned policy representations of 100 Monte Carlo runs. Figure S-2c illustrates an example state partitioning of one such
experiment, using the inferred local controllers depicted in
Fig. S-2e. The result reveals that all expert motion patterns
could be identified by our model. Note, however, that the
two figures correspond to a single Gibbs sample of the
process, which is not representative for the whole posterior
distribution. Averaging over several experiments as done in
Fig. S-2a and Fig. S-2b is not possible at the sample level due
to the varying dimensionality of the corresponding policy
representations (i.e. the number of learned controllers). Even
taking averages over samples of equal dimensionality is
not meaningful due to the multimodality of the posterior
distribution, which arises from the inherent symmetry of
the representation (i.e. interchanging two local controllers
together with their corresponding indices yields the same
model). Hence, averaging samples is possible only at the
prediction level.
Finally, Fig. S-2d depicts the posterior distribution of the
number of local controllers used by the model, which shows
a pronounced peak at the true number used by the expert.

C

C OMPUTATIONAL C OMPLEXITY

The overall computational cost of performing inference in
our model depends largely on two factors: the complexity
per Gibbs iteration and the mixing speed of the underlying
Markov chain. Each Gibbs iteration consists of up to three
stages: 1) sampling T categorical action variables {at } from
the set {1, . . . , |A|}, where T is the size of the demonstration set; 2) ddCRP model: sampling NS categorical state
assignments {ci } from the set {1, . . . , NS }, where NS is
e ); remaining models:
the number of states (i.e. |S| or |S|
sampling NS categorical partition assignments {zi } from
the set {1, . . . , K}, where K is the number of local controllers; 3) for non-collapsed models: sampling K Dirichletdistributed control parameters {θ k } on the (|A| − 1)simplex.
Collapsing the control parameters generally improves
the mixing speed of the chain (see Fig. 5 in the paper) but
requires that action variables belonging to the same cluster
be updated sequentially; hence, a non-collapsed strategy can

be advantageous for larger data sets. Sampling the variables
{at }, {θ k } and {zi } is computationally cheap because the
involved action likelihoods {T (s0 | s, a)} as well as the
neighborhood structure N (Potts model) and the similarity
values {f (di,j )} can be pre-computed. The most demanding
operation is the update of {ci }, which requires tracking
the connected components of the underlying ddCRP graph.
Using an appropriate graph representation, this can be done
in polylogarithmic worst case time [3].

R EFERENCES
[1] D. M. Blei and P. I. Frazier, “Distance dependent Chinese restaurant
processes,” The Journal of Machine Learning Research, vol. 12, pp.
2461–2488, 2011.
[2] C. E. Rasmussen, “Gaussian processes for machine learning.” MIT
Press, 2006.
[3] B. M. Kapron, V. King, and B. Mountjoy, “Dynamic graph connectivity in polylogarithmic worst case time,” in Proc. 24th Annual
ACM-SIAM Symposium on Discrete Algorithms, 2013, pp. 1131–1142.

