Recovering Markov Models from Closed-Loop Data

TL;DR

This paper develops an EM-based method to recover unbiased Markov models from data influenced by recommender systems that create feedback loops, addressing the challenge of hidden state dependencies.

Contribution

It introduces an EM algorithm to estimate transition matrices in Markov chains affected by feedback from recommender systems, accounting for the feedback loop.

Findings

The proposed method accurately estimates transition matrices in simulated experiments.

The approach effectively handles the feedback loop in Markov models.

Experimental results demonstrate the algorithm's robustness and applicability.

Abstract

Situations in which recommender systems are used to augument decision making are becoming prevalent in many application domains. Almost always, these prediction tools (recommenders) are created with a view to affecting behavioural change. Clearly, successful applications actuating behavioural change, affect the original model underpinning the predictor, leading to an inconsistency. This feedback loop is often not considered in standard so-called Big Data learning techniques which rely upon machine learning/statistical learning machinery. The objective of this paper is to develop tools that recover unbiased user models in the presence of recommenders. More specifically, we assume that we observe a time series which is a trajectory of a Markov chain ${R}$ modulated by another Markov chain ${S}$, i.e. the transition matrix of ${R}$ is unknown and depends on the current state of ${S}$. The transition matrix of the latter is also unknown. In other words, at each time instant, ${S}$ selects a transition matrix for ${R}$ within a given set which consists of known and unknown matrices. The state of ${S}$, in turn, depends on the current state of ${R}$ thus introducing a feedback loop. We propose an Expectation-Maximization (EM) type algorithm, which estimates the transition matrices of ${S}$ and ${R}$. Experimental results are given to demonstrate the efficacy of the approach.

Equations

$$P(R_t=j\mid R_{t-1}=i,S_{t-1}=s)=a_{ij}^R(s)$$

$$\begin{gathered} P(Q_t=(s^{\prime },r^{\prime })\mid Q_{t-1}=(s,r))= \\ P(R_t=r^{\prime }\mid R_{t-1}=r,S_{t-1}=s,\cancel{S_t=s^{\prime }})\cdot P(S_t=s^{\prime }\mid \cancel{R_{t-1}=r,}{S}_{t-1}=s)=a_{r{r}^{\prime }}^R(s)a_{s{s}^{\prime }}^S, \end{gathered}$$

$$P(S_t=j\mid S_{t-1}=i,R_{t-1}=r)=a_{ij}^S(\gamma (r)).$$

$$\begin{gathered} P(Q_t=(s^{\prime },r^{\prime })\mid Q_{t-1}=(s,r))= \\ P(R_t=r^{\prime }\mid R_{t-1}=r,S_{t-1}=s)\cdot P(S_t=s^{\prime }\mid R_{t-1}=r,S_{t-1}=s) \\ =a_{r{r}^{\prime }}^R(s)a_{s{s}^{\prime }}^S(\gamma (r)). \end{gathered}$$

$$\begin{split}\pi&=\pi^{S}\otimes\pi^{R},\qquad B=1_{S}\otimes I_{R},\\[4.30554pt] w_{ij}&=a^{S}_{\left\lceil i/R\right\rceil\,\left\lceil j/R\right\rceil}\bigl{(}\gamma(j_{R}(i))\bigr{)}\,a^{R}_{j_{R}(i)\,j_{R}(j)}\bigl{(}\left\lceil i/R\right\rceil\bigr{)},\end{split}$$

$$\begin{array}{rl}P(s_0,& \dots ,s_T,r_0,\dots ,r_T\mid \mu )\\ & ={\pi }_{r_0}^R{\pi }_{s_0}^S\prod _{t=1}^{T}P(s_t,r_t\mid s_{t-1},r_{t-1},\mu )={\pi }_{r_0}^R{\pi }_{s_0}^S\prod _{t=1}^T{a}_{s_{t-1}{s}_t}^S(\gamma (r_{t-1}))a_{r_{t-1}{r}_t}^R(s_{t-1})\end{array}$$

$$\begin{array}{rl}P& (r_0{r}_1\cdots r_T\mid \mu )\\ & =\sum _{s_0\in [S]\cdots s_T\in [S]}P(s_0,\dots ,s_T,r_0,\dots ,r_T\mid \mu )\\ & ={\pi }_{r_0}^R\sum _{s_0\cdots s_T}{\pi }_{s_0}^S\prod _{t=1}^T{a}_{s_{t-1}{s}_t}^S(\gamma (r_{t-1}))a_{r_{t-1}{r}_t}^R(s_{t-1})=:\ell (\mu )\end{array}$$

$${\alpha }_t(i):=P(S_t=i,R_0=r_0,\dots ,R_t=r_t\mid \mu )$$

$${\alpha }_t(j)=\sum _{i=1}^S{\alpha }_{t-1}(i)a_{ij}^S(\gamma (r_{t-1}))a_{r_{t-1}{r}_t}^R(i),j=1,\dots S,$$

$${\alpha }_t={\left(A^S(\gamma (r_{t-1}))\right)}^{⊺}\left(a_{r_{t-1}{r}_t}^R(:)\circ {\alpha }_{t-1}\right),$$

$$\ell(\mu)\!=\!\sum_{i=1}^{S}\!P(S_{T}=i,r_{0},\dotsc,r_{T}\mid\mu)\!=\!\sum_{i}\alpha_{T}(i)=\scalebox{1.05}{$\mathbf{1}$}^{\intercal}\alpha_{T}$$

$${\beta }_t(i):=P(r_{t+1},\dots ,r_T\mid S_t=i,r_t),$$

$${\beta }_{t-1}(j)=\sum _i{\beta }_t(i)a_{r_{t-1}{r}_t}^R(j)a_{ji}(\gamma (r_{t-1})),$$

$${\beta }_{t-1}=\left(A^S(\gamma (r_{t-1})){\beta }_t\right)\circ a_{r_{t-1}{r}_t}^R(:).$$

$$\alpha\mu+(1-\alpha)\nu=\Bigl{(}\alpha\pi^{R}+(1-\alpha)\rho^{R},\alpha\pi^{S}+(1-\alpha)\rho^{S},\\ \alpha A^{R}+(1-\alpha)B^{R},\alpha A^{S}+(1-\alpha)B^{S};\Gamma\Bigr{)}.$$

$$Q(\mu ,{\mu }^{\prime }):=\sum _{s_0\cdots s_T}P(s_0,\dots ,s_T,r_0,\dots ,r_T\mid \mu )\cdot \log P(s_0,\dots ,s_T,r_0,\dots ,r_T\mid {\mu }^{\prime }),$$

$$Q(\mu ,{\mu }^{\prime })=\log {\pi }_{r_0}^{R^{\prime }}\ell (\mu )+\sum _{i=1}^S\log {\pi }_i^{S^{\prime }}\sum _{j=1}^S{\xi }_1(i,j)+\sum _{i,j=1}^S\sum _{t=1}^T{L}_{ij}(r_{t-1},r_t){\xi }_t(i,j)$$

$${\xi }_t(i,j):=P(S_{t-1}=i,S_t=j,r_0,\dots ,r_T\mid \mu )$$

$${\xi }_t(i,j)={\alpha }_{t-1}(i)a_{r_{t-1}{r}_t}^R(i)a_{ij}^S(\gamma (r_{t-1})){\beta }_t(j),$$

$$\ell(\mu)\bigl{(}\log\ell(\mu^{\prime})-\log\ell(\mu)\bigr{)}\geq{Q(\mu,\mu^{\prime})-Q(\mu,\mu)}\,.$$

$${\pi }_i^{S^{\prime }}$$

$$a_{ij}^{S^{\prime }}(l)$$

$$a_{mn}^{R^{\prime }}(i)$$

$$\Vert f-g{\Vert }_{TV}=\underset{A\subseteq \Omega }{\max }f(A)-g(A)=\frac{1}{2}\sum _{\omega \in \Omega }|f(\omega )-g(\omega )|.$$

$$\Vert A^{R,\text{est}}(1)-A^{R,\text{true}}(1){\Vert }_{\text{stat}}:=\Vert {\psi }^{\text{est}}-{\psi }^{\text{true}}{\Vert }_{TV}.$$

$$\Vert A^{R,\text{est}}(1)-A^{R,\text{true}}(1){\Vert }_{\text{exp}}:=\sum _i{\psi }_i^{\text{true}}\Vert A_{i:}^{R,\text{est}}(1)-A_{i:}^{R,\text{true}}(1){\Vert }_{TV},$$

$$c_t$$

$$\begin{split}\hat{\alpha}^{\prime}_{t}&=\left(A^{S}(\gamma(r_{t-1}))\right)^{\intercal}\left(a^{R}_{r_{t-1}r_{t}}(:)\circ\hat{\alpha}_{t-1}\right)\\ \hat{\alpha}_{t}&=\frac{\hat{\alpha}^{\prime}_{t}}{\scalebox{1.05}{$\mathbf{1}$}^{\intercal}\hat{\alpha}^{\prime}_{t}}.\end{split}$$

$${\hat{\beta }}_{T-t}:={\beta }_{T-t}/c_t.$$

$$\begin{split}\hat{\xi}^{\prime}_{t}(i,j)&=\hat{\alpha}_{t-1}(i)a^{R}_{r_{t-1}r_{t}}(i)a^{S}_{ij}(\gamma(r_{t-1}))\hat{\beta}_{t}(j)\\ \hat{\xi}_{t}&=\frac{\hat{\xi}^{\prime}_{t}}{\scalebox{1.05}{$\mathbf{1}$}^{\intercal}\hat{\xi}^{\prime}_{t}\scalebox{1.05}{$\mathbf{1}$}}=\frac{{\xi}_{t}}{\scalebox{1.05}{$\mathbf{1}$}^{\intercal}\alpha_{T}}\end{split}$$

$${\pi }_i^{S^{\prime }}$$

Full text

1 Introduction

Our starting point for this paper is a frequently encountered problem that arises in the Smart Cities domain. Many decision support/recommender systems that are designed to solve Smart City problems are data-driven: that is data, sometimes in real time, is used to build models to drive the design of recommender systems. Almost always, these datasets are treated as if they were obtained in an open-loop setting, i.e. without recommender influence. However, this is rarely the case and frequently the effects of recommenders are inherent in datasets used for model building [18, 29, 8, 5]. This creates new challenges for the design of decision support/recommender systems under feedback. In particular, as engineers, we must take into account the fact that when we make a prediction, then this prediction affects the behaviour of operators [6] and this, in turn, changes the data set upon which the original model was built. Clearly, the aforementioned effect is related to classical closed-loop identification, which is itself a mature topic in both control and economics [31, 32, 12].

Notwithstanding this fact, and even though avoiding closed-loop effects in the design of recommender systems has been the subject of study in several Smart City applications [27, 26], the development of algorithms to identify models in closed loop remains a challenging problem in the context of Smart Cities. This is due to the fact that closed-loop questions that arise in Smart Cities are, for the most part, qualitatively different to those arising in other areas. For example, in control theory closed-loop questions typically arise in the context of deterministic and parametric models subject to noise, whereas in Smart Cities, typical problems are characterised by large-scale data sets that are generated by largely unknown stochastic processes.

Our objective is to consider one such problem class that arises in Smart City related research, where we seek to identify a user model, based on observations obtained when the user is acting under the influence of a recommender. A particular instance of such a problem arises in the automotive domain where drivers are characterised using Markovian models [15, 11], but where observations are obtained under the influence of recommenders acting on the driver. Motivated by such applications, we seek to develop methods to account for the effect of these recommender systems in data sets. More formally, we shall consider systems with the following structure: the process, which generates the data; the model, which represents the behaviour of the process; and the decision support tool, which intermittently influences the process. In our setup, data from the process is used to build the model. Typically, the model is used to construct a decision support tool which itself then influences the process directly. This creates a feedback loop in which the process, decision support tool and the model are interconnected in a complicated manner. As a result, the effect of the decision support tool is to bias the data being generated by the process, and consequently to bias any model that is constructed naively from the data.

To provide a little more context, and to return to the automotive example, we now illustrate such effects by means of the following application that we have developed in the context of our automotive research222https://www.youtube.com/watch?v=KUKxZZByIUM. Consider a driver who drives a car regularly. In order to design a recommender system for this driver we would like to build a model of his/her behaviour. For example, in order to warn the driver of, say, roadworks, along a likely route, we might use this model to predict the route of the driver.

A schematic of the proposed in-car architecture is depicted in Figure 1. The recommender uses a model of driver behaviour to issue intermittent recommendations. Observations of driver behaviour are then used to build a refined driver model which in turn is used as an input to the recommender system. Clearly, the effect of the recommender is to bias the driver model over time, thus eventually rendering the latter ineffective as an input to the recommender. The problems are exacerbated in many practical systems due to the presence of several unknown third-party recommender systems (Google Maps, Siri etc.), and by the fact that the driver model may operate from birth-to-death333 by that we mean that the driver always operates under the potential influence of a recommender. Thus, given any observation, we do not know whether the recommender is acting, or the driver. in closed loop. This latter fact makes it difficult, or impossible, to even estimate an initial model of driver behaviour. Clearly, in such applications it is absolutely necessary to develop techniques that extract the behaviour of the driver while under the influence of the feedback from a number of recommender systems. The results presented in this paper represent our first small step in this direction.

1.1 General Comments on Related Research Directions

Dealing with bias arising from closed-loop behaviour is a problem that has arisen in several application domains. In fact, in control theory, the related topic of closed-loop identification is considered to be a very mature area [22, 30]. Roughly speaking, this topic is concerned with building models of dynamic systems while they are being regulated by a controller. A related scenario arises in some adaptive systems when the controller itself is being adjusted on the basis of the dynamical systems model. As in our example, the controller action will bias the estimation of the model parameters. While many established techniques in control theory exist for dealing with such effects, these typically exploit known properties of the process noise and an assumed model structure to un-bias the estimates. Typically, structures such as ARMAX models are assumed to capture the nature of the system dynamics. Recent work on intermittent feedback [20] is also related in spirit to these approaches where control design techniques to deal with feedback loops that are intermittently broken are developed. Before proceeding it is worth noting that closed-loop effects have also been explored in the economics literature [33, 1, 31, 23].

More recently, several authors in the context of Smart Cities and recommender systems [8, 29], have realised that closed-loop effects represent a fundamental challenge in the design of recommender systems. In [8] the authors discuss the inherent closed-loop nature of data-sets in cities, and in [29] explicitly discuss the influence of feedback on the fidelity of recommender systems. As an example of a specific result, [29] presents an empirical technique for collaborative filtering to recover user rankings in the presence of a recommender under an assumed interaction model between user and recommender, which is similar in spirit to the aforementioned problem of constructing models from data possibly biased by a recommender, but does not consider sequential models.

The work reported here is closely aligned with stochastic models, unlike most of the approaches outlined in the first paragraph of Section 1.1. Specifically, in this paper we are interested in reconstructing Markov models that are operating under the influence of a recommender. To this end we assume that (i) recommenders and users can both be modelled by Markov chains, and (ii) recommendations are either accepted fully or have no influence at all, i.e. every decision is made by either the user or the recommender, never by a combination of the two. We note here that the Markovian assumption of user and recommender behaviour is convenient for many applications: for example, in the automotive domains [15, 17].

In this context, the present work is also related to classes of mixture and latent variable models, such as well-known hidden Markov (HMM) and mixture-of-experts (ME) models [21, 14]. In the latter case, a latent Markov chain selects from a set of parametrisations of a visible process, however there is no closed-loop modulation (i.e.  the modulated visible process is not allowed to in turn modulate the modulating latent process), and the visible process is static, subject to noise, not a Markov chain itself. Our work is most related to [9] but again with the important distinction of closed-loop modulation. A similar concept of regime switching time series models is used in econometrics [16]: these models allow parameters of the conditional mean and variance to vary according to some finite-valued stochastic process with states or regimes. However, the observations are assumed to be generated by a deterministic process with random noise, and the latent (switching) process is either a Markov chain independent of the past observations or is a deterministic function of the past observations. In contrast, we introduce the closed-loop modulation as discussed above. Yet another related model is Markov jump linear systems, see e.g. [7], where a (latent, autonomous) Markov chain selects the parameters of a (visible) dynamical system, whereas in our case, the visible part is a stochastic process on a discrete state space and it can modulate the latent process.

More specific technical comments to place our work in the context of reconstructing Markov models from data, and a brief discussion of practical issues including identifiability and convergence speed in terms of the number of samples are given below in Section 2.3 after the formal description of the proposed closed-loop Markov-modulated Markov chain models.

1.2 Preliminaries

Notation. To compactly represent discrete state spaces we write $[N]\mathrel{\mathop{:}}=\{1,\dotsc,N\}$. For a function $M:[N]\rightarrow\mathbb{R}^{m\times n}$ mapping such a discrete finite set to a set of matrices, we refer to each value $M(k)\in\mathbb{R}^{m\times n}$ as a page of $M$. Matrices will be denoted by capital letters, their elements by the same letter in lower case, and we denote the set of $n\times m$ row-stochastic matrices, i.e. matrices with non-negative entries such that every row sums up to 1, by $\mathcal{M}^{n\times m}$, and $\mathcal{M}^{m}\mathrel{\mathop{:}}=\mathcal{M}^{m\times m}$. For compatible matrices, $M\otimes N$ is the Kronecker product and $M\circ N$ denotes the Hadamard (or element-wise) product. A partition $\Gamma$ of $[N]$ is a set $\{\Gamma_{1},\dotsc,\Gamma_{p}\}$ such that $\Gamma_{i}\subseteq[N],\;\cup_{i}\Gamma=[N],\;\Gamma_{i}\cap\Gamma_{j}=\{\}$ $\forall i\neq j$. Each partition then also defines a membership function $\gamma:[N]\rightarrow[p]$ by $\gamma(i)\mathrel{\mathop{:}}=k\text{ such that }i\in\Gamma_{k}$. We write $P(W=w)$ for the probability of the event that a realisation of the discrete random variable $W$ equals $w$, and $P(W=w\mid V=v)$ the probability of that same event $W=w$ conditioned on the event $V=v$. We shall denote random variables by capital letters and, where appropriate, their realisations by the same letter in lower case. For convenience we will sometimes write $P(w\mid v)$ instead of $P(W=w\mid V=v)$ if there is no risk of ambiguity, and, for a set of parameters $\mu$ parametrising a probability distribution, $P(W=w\mid\mu)$ is taken to denote the probability of the event $W=w$ if the parameters are set to $\mu$.

Markov chains herein are sequences of random variables $\{X_{t}\}_{t}$ indexed by the time $t\in\mathbb{N}=\{0,1,2,\dotsc\}$. The realisation $x_{t}\in[N]$ of $X_{t}$ is the state of the Markov chain at time $t$, and $[N]$ is its state space. The probability distribution of $X_{0}$ is denoted by $\pi_{0}$, and the probability distribution of each following state is given by $P(X_{t}=j\mid X_{t-1}=i)=a_{ij}$; the matrix $A\in\mathcal{M}^{N}$ with entries $a_{ij}$ is the transition probability matrix.

2 Problem Statement and Model

As noted above, we assume that the driver and the recommender are Markovian. Given a possibly incomplete description of Markov chains modelling the recommender systems, and no knowledge of when these systems are engaged, our aim is to estimate the probability transition matrix of the Markov chain representing the driver, and the levels of engagement of each recommender, using only observed data. In what follows we formalize this setup, and give an expectation-maximization (EM) algorithm to estimate the parameters of the unknown driver model.

2.1 “Open-Loop” Markov-Modulated Markov Chains

Consider a Markov chain $\bm{R}$ with state space $[R]$ and state $r_{t}\in[R]$, in which the transition probabilities

[TABLE]

depend on a latent random variable $S_{t}$. We can say that the Markov chain is modulated by the random variable $S_{t}$, and if $S_{t}$ is itself the state of another Markov chain $\bm{S}$ with transition matrix $A^{S}$ and state space $[S]$, then we are dealing with a Markov-modulated Markov chain; Markov modulation is an established model in the literature on inhomogeneous stochastic processes, see e.g. [9].

Formally, the Markov-modulated Markov chain is defined by the tuple $\mu=(\pi^{R},\pi^{S},\allowbreak A^{R}(\cdot),A^{S})$, where $A^{R}:[S]\rightarrow\mathcal{M}^{R}$ and $\pi^{S},\pi^{R}$ denote the distributions of $S_{0}$ and $R_{0}$, respectively. That means for instance that if $\bm{S}$ only has a single state $s$, $\mu$ is a regular Markov chain with transition matrix $A^{R}(s)$ and initial probability $\pi^{R}$. We assume that we observe the state of $\bm{R}$, but not the state of $\bm{S}$.

Because the transition probabilities in the latent Markov chain $\bm{S}$ do not depend on the state of the visible chain $\bm{R}$, we refer to $\mu$ as an open-loop Markov-modulated Markov chain (ol3MC) to distinguish it from what follows. This models the case when the switching between the transition matrices $A^{R}$ occurs independently of the current state $r_{t}$ of $\bm{R}$.

The joint process ${Q_{t}}={(S_{t},R_{t})}$ has transition probabilities

[TABLE]

where the first cancellation means that the decision at time $t-1$ is not influenced by the state of the modulating random variable at time $t$, and the second cancellation follows from the open-loop assumption, i.e. that the modulating Markov chain $\bm{S}$ evolves independently of $\bm{R}$. The estimation of $\bm{S}$ and $\bm{R}$ for the case of continuous time ol3MCs has been discussed in [9].

Remark: We are dealing with the case when the data consists of a finite time series of observations of a single trajectory $(r_{0}r_{1}\dotsm r_{T})$ of the Markov chain $\bm{R}$ and no (estimate of the) distribution of $R_{t}$ is available. While if the distributions are available, standard methods of state-space identification apply, here the estimation of the parameters of $\mu$ requires statistical methods such as maximum likelihood estimation.

2.2 Closed-Loop Markov-Modulated Markov Chains

As a generalisation, we consider the case where $\bm{S}$ is dependent on the state of $\bm{R}$: that is, the probabilities of transitioning from one state $s$ to another state $s^{\prime}$ then do depend on what the current state $r_{t}$ is. We will be referring to this as a closed-loop Markov-modulated Markov chain or cl3MC for short. A cl3MC can be used to model that one transition matrix $A^{R}(s)$ might be more likely to be switched to in some regions of the visible state space $[R]$, or that switching can only occur when the system is in specific configurations. This is exactly the situation which arises in our automotive example, see Section 4.2.

Formally, we now also allow for the latent Markov chain $\bm{S}$ to be modulated by the current state of the visible chain $\bm{R}$. To keep the developments general, assume that – instead of one page of $A^{S}$ corresponding to *each *state of $\bm{R}$ – there is a partition $\Gamma=\{\Gamma_{1},\dotsc,\Gamma_{p}\}$ of $[R]$ such that there is a page in $A^{S}$ for each $\Gamma_{i}$. Hence, we now have $A^{S}:[p]\rightarrow\mathcal{M}^{S}$, with

[TABLE]

The open-loop case then corresponds to $\Gamma=\{[R]\}$ (i.e. $p=1$ and $\gamma(r)\equiv 1$) and the joint process ${Q_{t}}={(S_{t},R_{t})}$ has transition probabilities (compare to the open-loop formula above):

[TABLE]

Such a cl3MC is represented by a tuple $\mu=(\pi^{R},\allowbreak\pi^{S},\allowbreak A^{R}(\cdot),\allowbreak A^{S}(\cdot);\Gamma)$, where now, $A^{S}(\cdot)$ has pages, too.

To further illustrate the operation of a cl3MC model, Algorithm 1 details how a realization of the stochastic process described by it, i.e. a trajectory, is generated.

2.3 Relationship with Hidden Markov Models

There is a close relationship between closed-loop Markov modulated Markov chains and Hidden Markov Models (HMMs). Formally:

Proposition 1

$\mu=(\pi^{R},\pi^{S},A^{R}(\cdot),A^{S}(\cdot);\Gamma)$* defines the same visible process $\{R_{t}\}$ as the Hidden Markov Model $\lambda=(\pi,W,B)$, with*

[TABLE]

where $i,j=1,\dotsc,RS$, hence $W\in\mathcal{M}^{RS,RS}$.

Here, $j_{R}(k)\mathrel{\mathop{:}}=(k-1)(\bmod R)+1$ and $\left\lceil p\right\rceil\mathrel{\mathop{:}}=\inf\{k\in\mathbb{Z}\,\mid\,k\geq p\}$. The proof of Proposition 1 is mainly algebraic and can be found in [10].

Remark: While Proposition 1 maps a given cl3MC to an HMM which from the outside looks the same as, this mapping is not reversible: not every HMM represents a cl3MC, and most importantly, parameter estimation algorithms such as the standard Baum-Welch algorithm can not be used to estimate the parameters of a cl3MC, because they do not “respect the structure” of the matrix $W$: the HMM $\lambda$ is defined by $(RS)^{2}-RS+(RS-1)$ free parameters (the entries of $W$ and the entries of $\pi$ with the stochasticity constraints taken into account), whereas the corresponding cl3MC requires only $(pS-1)(S-1)+(SR-1)(R-1)$ parameters444 Note that $(pS-1)(S-1)+(SR-1)(R-1)<(RS)^{2}-1$ for $R+S>2$.. Hence, it is not possible to estimate the parameters of $\lambda$ and then compute the ones of $\mu$; instead, we develop an EM-algorithm to estimate the parameters of $\mu$ directly in Section 3.3.

**Identifiability: ** Given the close relationship between cl3MCs and HMMs outlined above, one should expect that identifiability issues for cl3MCs bear close resemblance to those of HMMs. By identifiability we mean the following: assume that $(r_{1}\dots r_{T})$ has been generated by the “true model” $\mu^{\mathrm{true}}$; under which conditions and in what sense will the estimate $\mu^{\mathrm{est}}$ converge to $\mu^{\mathrm{true}}$ if $T\to\infty$? For HMMs this question was partially answered in [24], namely it was shown that there is an open, full-measure subset $U$ of all HMMs, such that the sequence of estimates of the BW algorithm converges to $\lambda^{\mathrm{true}}$ (or a trivial permutation of it), provided the starting model is chosen within $U$, $\lambda^{\mathrm{true}}_{i}>\delta>0$ and $T\to\infty$. However, the structure of $U$ and convergence speed in terms of the number of samples were not described, and, to the best of our knowledge, these questions are still open.

For cl3MCs, similarly and trivially, any permutation of $A^{S}$ and the corresponding pages of $A^{R}$, which amounts to relabelling the hidden states $s$, yields the same visible process. However, there are examples of sets of HMMs $\lambda$, which are not permutations of each other, yet generate the same observable process; see [4, 13]. Interestingly, those examples involve the special case of partially observable Markov chains, a subclass of HMMs with emissions matrices $B$ having entries that are either $1$ or [math]. Comparing to Proposition 1, a cl3MC has close correspondence to an HMM of this class. This suggests that, in practice, the set of maximisers of the likelihood may be wider than the aforementioned set of permutations of $\mu^{\mathrm{true}}$; our numerical experiments in Section 4 also suggest that, in general, we cannot recover the true model $\mu^{\mathrm{true}}$, even up to trivial permutations, from observing only trajectories of $\bm{R}$. However, in the case of partial knowledge of elements of $A^{R}$ we can recover $A^{S}$ and the unknown portion of $A^{R}$. Hence, for the “driver-recommender” problem the proposed method is of practical value. Estimates of the minimum amount of prior knowledge necessary are the subject of future research.

3 Likelihood and Parameter Estimation

In this section we develop an iterative algorithm to estimate the parameters of a cl3MC $\mu=(\pi^{R},\allowbreak\pi^{S},\allowbreak A^{S}(\cdot),\allowbreak A^{R}(\cdot);\Gamma)$ given a sequence of observations $(r_{0}r_{1}r_{2}\dotsm r_{T})$, a partition $\Gamma=\{\Gamma_{1},\dotsc,\Gamma_{p}\}$ of $[R]$ and the size $S$ of the state space of $\bm{S}$. The derivation is close in spirit to the classical Baum-Welch (BW) algorithm (see e.g. [25] and the numerous references therein): our algorithm maximises at every iteration a lower bound on the likelihood improvement, and gives rise to re-estimation formulae (14) that utilise forward and backward variables which differ in subtle ways from the ones of the BW algorithm.

3.1 Likelihood of $\mu$, Forward- and Backward Variables

Since the estimate to be obtained is a maximum likelihood (ML) estimate, the efficient computation of the likelihood of a given cl3MC $\mu$ plays a central role in what follows. For a given $\mu$, the joint probability of sequences $(r_{0}r_{1}\dotsm r_{T})$ and $(s_{0}s_{1}\dots s_{T})$ being the trajectories of the visible Markov chain $\bm{R}$ and latent Markov chain $\bm{S}$ is

[TABLE]

where the last equality follows by (2). This allows us to compute the probability of observing a sequence $(r_{0}r_{1}\dotsm r_{T})$ given $\mu$ as follows:

[TABLE]

where $\mu\mapsto\ell(\mu)$ is the likelihood of the model $\mu$. Computation using this direct expression requires on the order of $2\times T\times S^{T}$ operations, and is hence not feasible for large $T$. Instead, we define the forward variable $\alpha_{t}$ with $S$ elements

[TABLE]

which can be computed iteratively as follows: $\alpha_{0}(j)=\pi^{S}_{j}\pi^{R}_{r_{0}}$ and

[TABLE]

or, in matrix form: $\alpha_{0}=\pi^{S}\pi^{R}_{r_{0}}$ and

[TABLE]

where the notation $a^{R}_{r_{t-1}r_{t}}(:)$ means a column vector of the $(r_{t-1},r_{t})$-elements of the matrix $A^{R}(k)$ as $k$ runs from $1$ to $S$.555Very much analogous to Matlab’s colon notation, or slicing in numpy. It follows that

[TABLE]

can be computed with on the order of $TS^{2}$ computations.

An analogous concept that will be required later is the backward variable

[TABLE]

which can also be computed via iteration: $\beta_{T}(j)=1$,

[TABLE]

or in matrix form: \beta_{T}=\scalebox{1.05}{\mathbf{1}} and

[TABLE]

3.2 Auxiliary Function $Q(\mu,\mu^{\prime})$

Let $\Lambda$ denote the set of all cl3MCs. $\Lambda$ is then bounded and convex if we define convex combinations of cl3MCs $\mu=(\pi^{R},\allowbreak\pi^{S},\allowbreak A^{R}(\cdot),\allowbreak A^{S}(\cdot);\Gamma)$ and $\nu=(\rho^{R},\allowbreak\rho^{S},\allowbreak B^{R}(\cdot),\allowbreak B^{S}(\cdot);\Gamma)$ as

[TABLE]

See [10] for details. Following [2], we define the auxiliary function $Q(\mu,\mu^{\prime})$ of $\mu,\mu^{\prime}\in\Lambda$ by

[TABLE]

where $s_{0},\dotsc,s_{T}$ run through the $S^{T}$ possible sequences of the latent state $S_{t}$. If parameters $\mu^{\prime}_{i}$ are zero where $\mu_{i}>0$, then we can have the case $P(s_{0},\dotsc,r_{T}\mid\mu^{\prime})=0$ and $P(s_{0},\dotsc,r_{T}\mid\mu)>0$; in this case $Q(\mu,\mu^{\prime})\mathrel{\mathop{:}}=-\infty$. If $P(s_{0},\dotsc,r_{T}\mid\mu)=P(s_{0},\dotsc,r_{T}\mid\mu^{\prime})=0$, we set $Q(\mu,\mu^{\prime})=0$ which amounts to setting $0\log(0)=0$.

The following lemma establishes a representation for $Q$ in terms of the elements of $\mu^{\prime}=(\pi^{R^{\prime}},\pi^{S^{\prime}},A^{S^{\prime}},A^{R^{\prime}};\Gamma)$:

Lemma 2

The function $Q(\mu,\mu^{\prime})$ can be rewritten as

[TABLE]

where $L_{ij}(m,n)\mathrel{\mathop{:}}=\log a^{S^{\prime}}_{ij}(\gamma({m}))+\log a^{R^{\prime}}_{mn}(i)$ and

[TABLE]

can be computed as follows:

[TABLE]

and the variables carrying a $\bullet^{\prime}$ constitute $\mu^{\prime}$.

The proof is relegated to B. Additionally, an application of Jensen’s inequality yields (see [10]):

Lemma 3

The improvement in log-likelihood satisfies the lower bound

[TABLE]

3.3 EM-Algorithm for Parameter Estimation

The algorithm proceeds by maximising the lower bound on the log-likelihood improvement set forth in (13) at every iteration.

It should be clear from (10) that the best estimate of $\pi^{R}$ is the $r_{0}$-th canonical Euclidian basis vector $e_{r_{0}}$. The remaining parameters of $\mu$ can be iteratively estimated by repeatedly applying the following theorem:

Theorem 1

The unique maximizer $\mu^{\prime}=M(\mu)$ of $Q(\mu,\cdot)$ is given by

[TABLE]

where $i,j=1,\dotsc,S$, $l=1,\dotsc,p$, and $m,n=1,\dotsc,R$.

The proof is given in B. Formulae (14) provide the basis for the EM-type parameter estimation algorithm for cl3MC $\mu$: in its $k$-th iteration, the E-step consists of computing $\xi_{t}(\cdot,\cdot)$ from the current estimate $\mu^{k}$, and the M-step yields an updated estimate $\mu^{k+1}=M(\mu^{k})$ with improved likelihood, see also the pseudocode in Algorithm 2 in the Appendix. Note that $M$ has a unique fixed point, $\mu^{\infty}=\lim_{k}\mu^{k}$ which is, at the same time, a stationary point (possibly a local maxima) of the likelihood, see [10] for details.

4 Examples

Here we illustrate the algorithm’s efficacy in two scenarios: first with synthetic data, i.e. data generated from a cl3MC, denoted $\mu^{\text{true}}$; and second, in a toy example of a practical application, estimation of driver behaviour. In both cases, we assume that one decision-maker, specifically the matrix $A^{R}(2)$, and $\Gamma$ are known. For implementation details, in particular how to avoid arithmetic underflow by scaling, and a pseudocode, see A and Algorithm 2 therein, and [10]; for details on the experimental procedures, see C.

4.1 Synthetic Data

To explore the relationship between estimation error and number of samples, repeated the following for several values $T\in[500,75000]$: $N_{e}=100$ cl3MCs with $R=20$, $S=2$ and $\Gamma=\{[R]\}$ (i.e. the open-loop case) were generated, and then a trajectory of length $T$ for each of them. We then ran the algorithm with random initial guesses $\pi^{S}$, $A^{S}$, and $A^{R}(1)$. The same was repeated for the same $N_{e}$ cl3MCs, only that now, $\Gamma$ was a randomly selected partition of order 2, so that $p=2$, and a second random page $A^{S}(2)$ was added to $A^{S}$. In both cases, we assume $A^{R}(2)$ and $\Gamma$ to be known. The modification to the algorithm is trivial: $A^{R}(2)$ is simply not re-estimated.

As illustrated in Figures 2 and 3, we recover $A^{S}$ and $A^{R}(1)$ to high accuracy for large enough $T$. “Accuracy” is hereby measured through statistical distances: since the transition matrices of Markov chains consist of probability distributions – row $i$ being the distribution of the state following $i$ – absolute or relative matrix norms are not a good measure of distance between Markov chains. Instead, we consider a statistical distance between the estimated and true probability distributions. One of the simplest such distances is the total variation (TV) distance (see e.g. [19, Ch. 4]), which is given by the maximal difference in probability for any event between two distributions. For probability distributions $f$ and $g$ over a discrete space $\Omega$, this is simply

[TABLE]

We consider here two applications of TV distance to Markov chains. The first is to take the TV distance between the stationary distributions, which concretely amounts to considering the subset $\rho$ of the state space $[R]$ such that $P(X_{t}\in\rho\mid A^{R,\text{est}}(1))-P(X_{t}\in\rho\mid A^{R,\text{true}}(1))$ is maximised (for large enough times $t$ such that the stationary distribution is reached). If we let $\psi^{\text{true}}$ and $\psi^{\text{est}}$ denote the stationary distributions, then

[TABLE]

However, this is a coarse measure: different Markov chains can have equal stationary distributions. Hence, the second metric incorporates the distance between the individual rows by considering the expectation (under the true stationary distribution $\psi^{\text{true}}$) of the TV distance between the estimated and the true row; this equals the sum of the distances between the true and estimated transition probabilities from all states $i$, weighted by the probability of being in state $i$:

[TABLE]

where $M_{i:}$ denotes the $i$-th row of matrix $M$.

The effect of $T$ on the accuracy is explored in Figure 2. The error appears to decay as a power of $T$, however this is simply an observation; a theoretical analysis of the sample complexity and decay rates is part of future work to be done.

For a representative value of $T=5\cdot 10^{4}$, Figure 3 drills down further into the experimental results; the distance for both introduced metrics is often below 10%, but we also observe severe outliers. Note that we show $A^{R}(1)$ only, the analysis and results for $A^{S}(\cdot)$ are analogous and are hence omitted.

4.2 A Model of Driver Behaviour

Recent research, e.g. [11, 17, 28, 15], suggests that Markov-based models are good approximations of driver behaviour and can be used e.g. for route prediction. Here, we illustrate how cl3MCs can be used to identify a driver’s preferences when some trips are planned by a recommender system, whose preferences are known, while the other trips are planned by the driver.

Specifically, consider the map in the left panel of Figure 4, which depicts a (very small toy) model of a driver’s possible routes from origin “O” to destination “D.” The houses, as an example, correspond to schools, that should be avoided in the hour before classes start and after classes end for the day, so there is a route past them and one around them. We assume that if a trip falls into that time frame, the recommender takes over and, with known probabilities, routes the driver either past or around each school; these probabilities make up $A^{R}(2)$. Otherwise, the driver follows his/her preferences, which constitute $A^{R}(1)$; this is the matrix we would like to estimate.

We generated $N_{e}=50$ sets of data by simulating $N_{t}=200$ trips on the graph shown in the right panel of Figure 4; this is the line graph of the map, where each road segment corresponds to a node, and an edge goes from node $i$ to node $j$ iff it is possible to turn into road segment $j$ from $i$. Each trip has a probability of $p_{r}=0.3$ to be planned by the recommender. If a trip was planned by the recommender (resp. driver), a trajectory was generated by a Markov chain with transition matrix $A^{R}(2)$ (resp. $A^{R}(1)$) originating in node 1 and terminating when returning to node 1.

For estimation in the cl3MC framework, all trips are then concatenated to form one long trajectory and $A^{S}$ and $A^{R}(1)$ are estimated for an ol3MC, i.e. for $\Gamma=[9]$. $A^{R,\text{est}}(1)$ is then an estimate of the driver preferences. The results are shown in the first column of Figure 5 and are satisfactory already; however, we can leverage the closed-loop framework to include the additional knowledge that the the decision maker (i.e. the page of $A^{R}$ used) can only change after a trip is finished. Because the decision which page of $A^{R}$ to use at time $t$ is made at $t-1$, see (1), this means we have to allow for the state of $\mathbf{S}$ to change on the road segments prior to reaching the destination. We hence let \Gamma=\bigl{\{}\{8,9\},\{1,\dotsc,7\}\bigr{\}} and $A^{S}(2)=I_{2}$. $A^{S}(1)$ needs to be identified. The results are shown in Figure 5.

Additionally, we can interpret $\psi_{2}^{S,\text{est}}$, the second element of the stationary distribution of $A^{S,\text{est}}(1)$ as an estimate of $p_{r}$. For the open-loop case, we obtain $0.352(\pm 0.0583)$, whereas the cl3MC estimation yields $0.295(\pm 0.044)$.

5 Concluding Remarks

We consider the identification of user models acting under the influence of one or more recommender systems. As we have already discussed, actuating behavioural change affects the original model underpinning the predictor, leading to an biased user models. Given this background, the specific contribution of this paper is to develop techniques in which unbiased estimates of user behaviour can be recovered in the case where recommenders, users, and switching between them can be parameterised in a Markovian manner, and where users and recommenders form part of a feedback system. Examples are given to present the efficacy of our approach.

Acknowledgements: The authors would like to thank Ming-Ming Liu and Yingqi Gu (University College Dublin) for help with the numerical experiments, and Giovanni Russo and Jakub Mareček (IBM Research Ireland) for valuable discussions.

This work has been conducted within the ENABLE-S3 project that has received funding from the ECSEL joint undertaking under grant agreement NO 692455. This joint undertaking receives support from the European Union’s HORIZON 2020 Research and Innovation programme and Austria, Denmark, Germany, Finland, Czech Republic, Italy, Spain, Portugal, Poland, Ireland, Belgium, France, Netherlands, United Kingdom, Slovakia, Norway.

Robert Shorten was also partially supported by SFI grant 16/IA/4610.

Appendix

Appendix A Scaling Issues

Since the computations of $\alpha_{t}$ and $\beta_{T-t}$ according to (6) and (8) involve multiplications of on the order of $2t$ numbers less than 1, for large $T$, they will be close to, or below, machine precision. The re-estimation (14) then requires division of very small numbers, which of course should be avoided. To mitigate these issues, scale $\alpha_{t}$ to sum up to 1:

[TABLE]

The update can be done in two compact steps:

[TABLE]

Also note that \log c_{t}=\log c_{t-1}+\log\scalebox{1.05}{\mathbf{1}}^{\intercal}\hat{\alpha}^{\prime}_{t} and \log\ell(\mu)=\log(\scalebox{1.05}{\mathbf{1}}^{\intercal}\alpha_{T})=\log(c_{T}); the likelihood $\ell(\mu)$ is not computed anymore, only the log-likelihood. Since the backwards variables $\beta_{T-t}$ can be expected to be of similar order as $\alpha_{t}$ they are scaled using the same scaling factors $c_{t}$:

[TABLE]

We then compute the scaled versions of $\xi_{t}(i,j)$ in similar fashion

[TABLE]

and note that upon substituting \xi_{t}=(\scalebox{1.05}{\mathbf{1}}^{\intercal}\alpha_{T})\hat{\xi}_{t} in (14), \scalebox{1.05}{\mathbf{1}}^{\intercal}\alpha_{T} cancels everywhere, and we arrive at the rescaled re-estimation equations

[TABLE]

For a more detailed derivation, see again [10].

Appendix B Proofs

Proof 1** (Lemma 2)**

Formula (12) follows directly from the definitions of $\alpha_{t-1}$, $a^{R}_{mn}$, $a^{S}_{ij}$ and $\beta_{t}$. Now let summation indices $s_{\tau}$ always run from 1 to $S$. First, we rewrite the second line of (9): given the definition of $L_{ij}(k,\ell)$ and (3) we get

[TABLE]

We substitute this back into (9) to get $Q(\mu,\mu^{\prime})=Q_{1}+Q_{2}+Q_{3}$, with

[TABLE]

for $Q_{2}$ and $Q_{3}$, note the “marginalisation”

[TABLE]

where the sum runs over the states of the latent chain $\bf S$ before time instant ${t-1}$ and after time instant $t$. Then

[TABLE]

and

[TABLE]

This completes the proof.

Proof 2** (Theorem 1)**

From the remark before the theorem, it should be clear that we can ignore the first term in (10) in the maximisation. Consider $\mu\in\Lambda$ as fixed, and define $\widetilde{\mu}\mapsto W(\widetilde{\mu}):=Q(\mu,\widetilde{\mu})$. We claim that $W$ has the unique global maximum point $\mu^{\prime}\in\Lambda$. Note that if $\mu\in\Lambda$ then it may have zero components, say $\mu_{i_{1}}=\dots=\mu_{i_{d}}=0$. Then the logarithms of the corresponding components $\tilde{\mu}_{i_{1}}\dots\tilde{\mu}_{i_{d}}$ of $\tilde{\mu}$ in $W$ are multiplied by [math] so that these components do not change $W$. However, if we fix all the components of $\tilde{\mu}$ but $\tilde{\mu}_{i_{1}}$ and any other component $\tilde{\mu}_{k}$ such that $k\not\in\{i_{1}\dots i_{d}\}$, and $k$ is such that $\tilde{\mu}_{i_{1}}$ and $\tilde{\mu}_{k}$ are in $\pi^{R}$, or are in the same row of $A^{R}$ or $A^{S}$, then increasing $\tilde{\mu}_{i_{1}}$ will decrease $\tilde{\mu}_{k}$ (to meet the stochasticity constraints). As a result, the $\log\tilde{\mu}_{k}$ will decrease causing $W$ to decrease. Hence, the maximum of $W$ is attained in the set $\widetilde{\Lambda}\mathrel{\mathop{:}}=\{\tilde{\mu}\in\Lambda:\tilde{\mu}_{i_{1}}=\dots=\tilde{\mu}_{i_{d}}=0\}$. Let $\widetilde{W}$ be the restriction of $W$ to $\widetilde{\Lambda}$. Now $\widetilde{W}$ is a conical sum of logarithms of all $R+S+SR^{2}+pS^{2}-d$ independent components of $\tilde{\mu}$, hence strictly concave function on a convex compact set $\widetilde{\Lambda}$. Hence, $\widetilde{W}$ has a unique maximum point in $\widetilde{\Lambda}$ which coincides with the unique global maximum point $\mu^{\prime}$ of $W$ in $\Lambda$.

Let us prove that $\mu^{\prime}=M(\mu)$. As noted above, $\mu_{i_{1}}=\dots=\mu_{i_{d}}=0$ implies that $\mu^{\prime}_{i_{1}}=\dots=\mu^{\prime}_{i_{d}}=0$, and we stress that the same property holds true for $M(\mu)$: If $\pi^{S}_{i}$ is [math] it follows from (6) that $\alpha_{0}(i)=0$ and from (12) we get $\xi_{1}(i,j)=0$. By (14a), $\pi^{S^{\prime}}_{i}=0$ as well. If we have $a^{S}_{ij}(l)=0$, then again from (12), we see that $\xi_{t}(i,j)=0$ for all $t$ with $r_{t-1}=l$, and (14b) yields $a_{ij}^{S^{\prime}}(l)=0$. Similarly, $a^{R}_{mn}(i)=0$ leads to $\xi_{t}(i,j)=0$ whenever $r_{t-1}=m$ and $r_{t}=n$, independently of $j$, so $a^{R^{\prime}}_{mn}(i)=0$, too. Hence $M(\widetilde{\Lambda})=\widetilde{\Lambda}$, and $\pi^{S^{\prime}}_{i}$, $a^{S^{\prime}}_{ij}(l)$ and $a^{R^{\prime}}_{mn}(i)$ defined by (14) are positive if the corresponding components of $\mu$ are. On the other hand, $\mu^{\prime}_{k}>0$ if $\mu_{k}>0$ as otherwise $W(\mu^{\prime})=-\infty$, and so the gradient of $\widetilde{W}$ is well-defined at $\mu^{\prime}$. In fact, for positive $\tilde{\pi}^{S}_{i}$, $\tilde{a}^{S}_{ij}(l)$ and $\tilde{a}^{R}_{mn}(i)$:

[TABLE]

By e.g. [3, p.113, Prop. 2.1.2], it is necessary and sufficient for $\mu^{\prime}$ to satisfy the inequality

[TABLE]

We stress that the r.h.s. of (21) is independent of $\mu^{\prime}$ as

[TABLE]

Since $\pi^{S^{\prime}}_{i}$, $a^{S^{\prime}}_{ij}(l)$ and $a^{R^{\prime}}_{mn}(i)$ defined by (14) are positive if the corresponding components of $\mu$ are so, the gradient of $\widetilde{W}$ is well-defined at $M(\mu)$. Take any $\widetilde{\mu}\in\widetilde{\Lambda}$ and compute:

[TABLE]

so that, by stochasticity constraint, we get [10]:

[TABLE]

Hence, $M(\mu)$ defined by (14) satisfies (21) with equality for any $\widetilde{\mu}\in\widetilde{\Lambda}$. This completes the proof.

Appendix C Experimental details

The detailed steps in generating the data in Section 4.1 are as follows. We consider ${T}\in\{500,\allowbreak 1000,\allowbreak 3000,\allowbreak 5000,\allowbreak 7500,\allowbreak 10000,\allowbreak 12500,\allowbreak 15000,\allowbreak 20000,\allowbreak 25000,\allowbreak 30000,\allowbreak 50000,75000\}$ and begin at $i=0$. Then, for each of the $N_{T}=13$ values of $T$, the following is repeated $N_{e}$ times:

1. generate a pair of $R\times R$ row-stochastic matrices $\left(A^{R}(1),A^{R}(2)\right)$ by selecting one “dominant element” per row and setting it to a random number uniformly distributed in $[0.5,1]$ (this is done to ensure that the two matrices in a pair are sufficiently different) and filling the remaining elements with uniformly random numbers;

2. generate a $S\times S$ row-stochastic matrix $A^{S}(1)$ with uniformly random entries;

3. generate initial probability vectors $\pi^{S}\in[0,1]^{S}$ and $\pi^{R}\in[0,1]^{R}$ with uniformly random entries;

4. generate a trajectory of length $T$ from the cl3MC $\mu_{o}=(\pi^{S},\pi^{R},A^{S}(1),A^{R};[R])$ (this is the open-loop case);

5. estimate parameters of $\mu_{o}$ by running Algorithm 2 (initialized with uniformly random guesses for unknown parameters) until convergence and record the distance measures defined in Section 4.1;

6. generate an additional $S\times S$ row-stochastic matrix $A^{S}(2)$ with uniformly random entries;

7. generate a random partition $\Gamma$ of $[R]$ by first randomly permuting $[R]$ and then splitting it after a random index between $1$ and $R-2$;

8. generate a trajectory of length $T$ from the cl3MC $\mu_{c}=(\pi^{S},\pi^{R},A^{S},A^{R};\Gamma)$ (this is the open-loop case);

9. estimate parameters of $\mu_{c}$ by running Algorithm 2 (initialized with uniformly random guesses for unknown parameters) until convergence and record the distance measures defined in Section 4.1.

Note that we assume $A^{R}(2)$ to be known and it is not estimated. Concretely that means that the initial guess of $A^{R}(2)$ is set to the true value of $A^{R}(2)$ and $A^{R}(2)$ is excluded from the reestimation steps in Algorithm 2. The parameters used in Algorithm 2 are $\varepsilon=10^{-5}$, a typical, empirical choice for relative tolerances, and $N_{max}=2000$, at which point the algorithm has typically long converged (in our experiments, $N_{max}$ was never reached). Once the experiment terminates, we have $2N_{e}N_{T}=2600$ (resp.  $3N_{e}N_{T}=3900$) different values for each of the distance metrics for $A^{R}$ (resp. $A^{S}$), which are visualized in Figures 2 and 3.

For Section 4.2, the matrices $\left(A^{R}(1),A^{R}(2)\right)$ are fixed to the ones indicated in Figure 4. One run then consisted of repeating $N_{t}$ times the following:

1. draw a number from the uniform distribution on $[0,1]$ and if it is less than $p$, set $s=1$, else $s=2$;

2. using the transition probabilities in $A^{R}(s)$, generate a trajectory from initial state $O$ until state $D$ is reached.

We then concatenated all of those trajectories, identifying states $O$ and $D$ as state 1, e.g. if there were only two trajectories, $(O,2,6,5,9,D)$ and $(O,7,4,5,8,D)$ they are concatenated to $(1,2,6,5,\allowbreak 9,1,7,4,5,8,1)$. Algorithm 2 is then run on that single trajectory with the same parameters as above, again assuming $A^{R}(2)$ as known, two times: once assuming open-loop, i.e. $\Gamma=[9]$, and once assuming closed-loop with \Gamma=\bigl{\{}\{8,9\},\{1,\dotsc,7\}\bigr{\}} and $A^{S}(2)=I_{2}$. Then the distance metrics between estimates and true values are computed. This concluded one run; we collected $N_{e}$ runs, yielding $2N_{e}=100$ numbers for each $A^{R}(1)$ and $A^{S}(1)$. The numbers for $A^{R}(1)$ are shown in Figure 5.