Indirect Observation of State and Transition Probabilities

Abstract

A wide range of systems exhibit stochastic transitions between different states that may be hidden from direct observation. Nevertheless, if the states are coupled to a signal, observation of the signal can provide necessary information to infer the state and switching characteristics. Here we explore a simple hidden Markov model with an observable Poissonian distributed count signal. Determining the parameters of this system from the signal can be difficult in the high-noise regime with non-Bayesian methods. However this system yields a simple Bayesian network description, and variable independencies allow the problem to be formulated in a way that allows tractable inference of the parameters just from the time series. This is an informative demonstration of Bayesian techniques, and in particular the interplay between modelling a system and the process of inference.

1. Introduction

One of the really satisfying aspects of Bayesian inference is the structured interplay between modelling and observable data. To begin an analysis, a clear model of the system and noise sources is required. The model informs how to infer parameters from correlations in the observed data. In turn, comparison with experiment might suggest extensions or modifications of the model and lead to more rounds of analysis. We demonstrate this interplay in this chapter with a simple hidden Markov model, that is a distillation of the analysis presented [1] which was applied to quantum emitters. The original system in question had just two states, labelled on and off, with a probabilistic transition between these states that does not depend on the past of the system. We extend the modelling to an arbitrary number of states and transition probabilities. In the scenario we consider, the system state and transitions are hidden from the observer. What is observable is a Poissonian distributed signal whose strength depends on the system state. The task is then to determine the hidden parameters from the observed signal, such as the system state over time, or the switching probabilities.

A variety of systems exhibit stochastic switching between states, sometimes called telegraph noise [2], and the first step to controlling the system might be to understand the characteristics of the switching. An example of such a system is a quantum emitter, where the switching behaviour is called blinking due to large intermittent changes in the observed fluorescence. Such emitters are valuable resources for engineered quantum systems, and the presence of blinking limits their applicability. If fact, most quantum light emitters exhibit intermittency in their fluorescence under continuous excitation, including diamond colour centres [3, 4, 5, 6, 7], quantum dots [8, 9, 10, 11, 12], nanowires [9, 13], nanorods [13], organic semiconductors [14], molecules [15, 16, 17], and other systems [18, 19]. Key insights can be distilled from the switching rates and how these depend on relevant parameters, and yet these rates are often the hardest to extract from the raw time series of photon counts in particular for low-intensity light signals or noisy data.

Blinking typically leads to step-like switches in the fluorescence time trace, as illustrated in Figure 1(a). The most common method used to analyse the on and off states from a blinking time series is threshold analysis [3, 20, 21], illustrated for simulated data in Figure 1(a). However, the choice of threshold intensity is essentially arbitrary, and it can significantly influence the statistics of the on and off states [22, 23, 24, 25]. The threshold technique becomes difficult to use for emitters fluorescing at low signal-to-noise ratios as depicted in Figure 1(b), where the distributions of counts from the on and off states to overlap. Figure 1(c) shows that generalising this situation to three (or more) possible states typically makes it less likely that suitable thresholds can be found to distinguish all of the states.

Figure 1.

Not all switching time traces suit analysis in terms of a threshold. (a) Simulated data with clear on and off states. The dotted line represents a threshold, which is used to distinguish the on and off states. (b) An example of the same model with the same switching rates where it would be difficult to do a threshold analysis. (c) A model with three states with no clear distinction.

By taking a Bayesian inference approach we can extract state switching probabilities or rates directly from the time-series of the observable parameter, without resorting to threshold analysis. In order to not over complicate the analysis we will focus on a regime where the state switching rates are slow in comparison to the detection intervals. The analysis can also be extended to the continuous switching regime as demonstrated in [1], either as a sequence of more and more sub-intervals or in the fully continuum limit.

In Section 2 we analyse a switching process where the system is assumed to be either in a number of possible states $\left\{a,b,c,\dots \right\}$ for the whole detection interval. In Section 3 we infer the switching rates from observed data. Finally in Section 4 we use the time trace to infer the system state itself.

2. Modelling

The model for the system has some number of hidden states, $\left\{a,b,c,\dots \right\}$, and it switches randomly from a state $x$ to a state $y$ with probability of $T_{\mathit{xy}}$. The system state is only indirectly observable by a signal, that is itself Poissonian distributed with a rate $R_s$ that depends on the system state $s$. The Markov chain for this model is depicted in Figure 2 for three states.

Figure 2.

A discrete time Markov process modelling the evolution of the system between three states. In this case $T_{\mathit{xy}}$ are the switching probabilities from state $x$ to state $y$. With three states there are potentially six different transition probabilities to track.

Such a model can be also used to model a background Poissonian noise process that affects all rates by simply adding the background rate to all other rates, as a combination of Poissonian processes is also Poissonian. In the physical example of quantum emitters, the background rate (known as dark-counts) is due to noise in the detector or electronics. Similarly, while the quantum emitter is in the on state and is fluorescing, it is being driven by a strong classical pump so the statistics of the photon counts are also well modelled by a Poissonian process. The observed counts are then Poissonian with a rate of either the background rate or the background rate and the fluorescing rate combined, both of which are Poissonian.

We’ll write all the state dependent rate parameters as a single variable for brevity in the probability statements $R=R_a{R}_b{R}_c\dots$. The state at the beginning of time-step $t$ (and throughout the time-step) is $s_{t-1}$, and the probability of seeing $c_t$ counts integrated over the detection window is given by a Poisson distribution with a rate determined by this state,

$$P\left(c_t|s_{t-1}\mathit{RI}\right)=c_t^{\lambda }exp\left(-\lambda \right)/c_t!$$E1

where $\lambda =R_{s_{t-1}}$, $s_{t-1}\in \left\{a,b,c,\dots \right\}$, and $I$ tags the background information for the model.

Clearly, if we knew the state of the emitter at each data point it would be a simple matter to infer the switching rates. Unfortunately these states are not directly observable, and worse still, because the states are not observed the switching probabilities become dependent on the entire history of the count data. This makes the inference considerably more involved.

We can summarise the dependencies amongst the variables by the Bayesian network (BN) [26] shown in Figure 3. The BN can be used to determine conditional independencies between the problem variables, using the property of $d$-separability. In the BN, $s_{t-1}$ represents the state of the emitter for the $t^{\mathrm{th}}$ detector interval. The key variable independencies are the following

Figure 3.

A partial Bayesian network representing the joint probability distribution of problem parameters. The nodes inside the central dark regions are the visible variables, those outside are hidden parameters. Each dark region represents a detection interval and the visible variables are the accumulated counts for the entire detection interval. We are estimating rate parameters $R_x\in R$, and transition probabilities $T_{\mathit{xy}}\in T$, for $x,y\in \left\{a,b,c,\dots \right\}$ given the observed counts $c_t$. Enough of the full network is drawn to be able to easily determine the variable independencies.

$$T\perp R\mid I$$E2

$$s_t\perp s_u\mid s_{t-1}\mathrm{\Omega }\left(u<t-1\right)$$E3

$$c_t\perp c_u\mid s_{t-1}\mathrm{\Omega }u\ne t$$E4

$$c_t\perp s_v\mid s_{t-1}\mathrm{\Omega }v\ne t-1$$E5

where $\mathrm{\Omega }=TRI$ for compactness. Note that ${\sum }_y{T}_{\mathit{xy}}=1$ for all $x$, that is all the exit probabilities should add up to one. We will make use of these independencies in the following section.

Generating sample data from the model is straight forward. The initial system state is set at some arbitrary state, say $s_0=a$. Then at each time step $t$ the next state is chosen according to the probabilities $\left\{T_{s_{t}x}\right\}$ for $x\in \left\{a,b,c,\dots \right\}$. This does a random walk between states following Figure 2. Each time step is “observed” by Possonian generated counts$c_t$ whose rate is determined by the state: $R_{s_t}$.

3. Inference

Given the observed counts, the quantity we want to infer is $P\left(T|\overrightarrow{c}I\right)$ where $\overrightarrow{c}\equiv c_1{c}_2\dots c_N$ is the set of count data obtained over $N$ detector intervals, and $T$ are all the state transition probabilities. Using Bayes’ rule, the posterior probability distribution of $T$ can be written as,

$$P\left(T|\overrightarrow{c}I\right)=\frac{P\left(T|I\right)P\left(\overrightarrow{c}|TI\right)}{P\left(\overrightarrow{c}|I\right)}$$E6

The rates of counts $R$ can be added by marginalisation,

$$P\left(T|\overrightarrow{c}I\right)=\int \mathit{dR}\frac{P\left(T|I\right)P\left(\overrightarrow{c}R|TI\right)}{P\left(\overrightarrow{c}|I\right)}$$E7

$$=\frac{1}{P\left(\overrightarrow{c}|I\right)}\int \mathit{dR}P\left(\mathit{TR}|I\right)P\left(\overrightarrow{c}|\mathrm{\Omega }\right).$$E8

Where the integral is over all the rates $\int \mathit{dR}=\int {\mathit{dR}}_a\int {\mathit{dR}}_b\int {\mathit{dR}}_c\dots$ for some suitable ranges. Without observing the counts, the independency in (2) implies we can factor $P\left(TR|I\right)=P\left(T|I\right)P\left(R|I\right)$. We will go further and take the probability of all parameters as constant over some initial range (i.e. no prior information) so that

$$P\left(T|\overrightarrow{c}I\right)=\frac{1}{\mathcal{N}}\int \mathit{dR}P\left(\overrightarrow{c}|\mathrm{\Omega }\right)$$E9

and determine the normalisation factor $\mathcal{N}$ at the end. Note that $\mathrm{\Omega }=TRI$ so these probabilities still depend on the rates and transition probabilities.

If all the states $\overrightarrow{s}\equiv s_0{s}_1\dots s_{N-1}$ were known together with the parameters $R$ and $T$, the probability of all the counts could again be easily determined. This suggests an approach to the problem by adding these parameters through marginalisation,

$$P\left(\overrightarrow{c}|\mathrm{\Omega }\right)=\sum _{\overrightarrow{s}}P\left(\overrightarrow{c}\overrightarrow{s}|\mathrm{\Omega }\right)=\sum _{\overrightarrow{s}}P\left(\overrightarrow{c}|\overrightarrow{s}\mathrm{\Omega }\right)P\left(\overrightarrow{s}|\mathrm{\Omega }\right)$$E10

where we abbreviate the summations as ${\sum }_{\overrightarrow{s}}={\sum }_{s_0}\dots {\sum }_{s_{N-1}}$. Using the independencies (4) and (5), $P\left(\overrightarrow{c}|\overrightarrow{s}\mathrm{\Omega }\right)$ can be simplified to

$$P\left(\overrightarrow{c}|\overrightarrow{s}\mathrm{\Omega }\right)=\prod _{t=1}^{N}P\left(c_t|s_{t-1}\mathrm{\Omega }\right),$$E11

and each term is determined by (1).

The remaining $P\left(\overrightarrow{s}|\mathrm{\Omega }\right)$ term cannot be simply factorised over the states despite being a Markov chain, as observing $R$ and $T$ introduces possible dependencies. We can however expand using the product rule and simplify by making use of the independencies in (3):

$$P\left(\overrightarrow{s}|\mathrm{\Omega }\right)=\prod _{t=1}^{N}P\left(s_t|s_{t-1}\mathrm{\Omega }\right)P\left(s_0|\mathrm{\Omega }\right),$$E12

where we have chosen to expand in temporal order.

Finally, the inference becomes

$$P\left(T|\overrightarrow{c}I\right)=\frac{1}{\mathcal{N}}\int \mathit{dR}\sum _{\overrightarrow{s}}\prod _{t=1}^{N}P\left(c_t{s}_t|s_{t-1}\mathrm{\Omega }\right)P\left(s_0|\mathrm{\Omega }\right)$$E13

where we have used (5) to write the joint distribution between $c_t$ and $s_t$.

The problem with (13) is that the sum over $\overrightarrow{s}$ contains $m^N$ terms if there are $m$ states, each of which has $N$ products. This will rapidly become intractable as the size of the data grows. Fortunately it’s possible to rewrite (13) as a single term with $N$$m\times m$ matrix products. Consider the following matrix written for a three-state system:

$$Q_t=\left[\begin{array}{ccc}P\left(c_t{s}_t=a|s_{t-1}=a\mathrm{\Omega }\right)& P\left(c_t{s}_t=a|s_{t-1}=b\mathrm{\Omega }\right)& P\left(c_t{s}_t=a|s_{t-1}=c\mathrm{\Omega }\right)\\ P\left(c_t{s}_t=b|s_{t-1}=a\mathrm{\Omega }\right)& P\left(c_t{s}_t=b|s_{t-1}=b\mathrm{\Omega }\right)& P\left(c_t{s}_t=b|s_{t-1}=c\mathrm{\Omega }\right)\\ P\left(c_t{s}_t=c|s_{t-1}=a\mathrm{\Omega }\right)& P\left(c_t{s}_t=c|s_{t-1}=b\mathrm{\Omega }\right)& P\left(c_t{s}_t=c|s_{t-1}=c\mathrm{\Omega }\right)\end{array}\right],$$E14

and vector

$${\overrightarrow{D}}_0=\left[\begin{array}{c}P\left(s_0=a|\mathrm{\Omega }\right)\\ P\left(s_0=b|\mathrm{\Omega }\right)\\ P\left(s_0=c|\mathrm{\Omega }\right)\end{array}\right]$$E15

then

$$\sum _{\overrightarrow{s}}\prod _{t=1}^{N}P\left(c_t{s}_t|s_{t-1}\mathrm{\Omega }\right)P\left(s_0|\mathrm{\Omega }\right)={\overrightarrow{1}}^T\prod _{t=1}^N{Q}_t{\overrightarrow{D}}_0$$E16

where the product on the right hand side is read as decreasing in $t$ to the right. The equivalence of (16) to (13) is readily verified by expanding a few terms out. The matrix multiplication will sum over the columns of $Q_t$ which is a summation over the states $s_{t-1}$. Each successive multiplication sums over another state and the final multiplication by the row vector ${\overrightarrow{1}}^T$ (for the three state system this is $\left[111\right]$) will sum over $s_N$. The generalisation of $Q_t$ and ${\overrightarrow{D}}_0$ to more than three states is straight forward. With this equivalence, the inference reads

$$P\left(T|\overrightarrow{c}I\right)=\frac{1}{\mathcal{N}}\int \mathit{dR}{\overrightarrow{1}}^T\prod _{t=1}^N{Q}_t{\overrightarrow{D}}_0.$$E17

There is one outstanding issue in calculating (17), and that is that care must be taken to avoid underflow or overflow. The normalisation factor includes $P\left(\overrightarrow{c}|I\right)$ which is the probability of the observed counts given no knowledge of the parameters or transition probabilities. While it could be expanded by adding parameters and states it would be difficult to calculate. Without it, the computation of (17) quickly leads to overflow. The key to controlling this normalisation issue is to observe that (17) is similar to a time evolving state

$$\prod _{t=1}^N{Q}_t{\overrightarrow{D}}_0=Q_N\dots Q_3{Q}_2{Q}_1{\overrightarrow{D}}_0=Q_N\dots Q_3{Q}_2{\overrightarrow{D}}_1=Q_N\dots Q_3{\overrightarrow{D}}_2.$$E18

So we can compute the “state” $D_t$ for a multi-dimensional grid of parameter values that spans the inference space, and at each time step re-normalise since the sum of the probabilities over the entire grid should be 1. In this way the normalisation can be controlled.

As a demonstration of the algorithm, Figure 4 shows the marginals of the posterior distributions for all the exit probabilities from each state inferred from the data of Figure 1(c). The count rates $R_a,R_b$ and $R_c$ where assumed known for the purpose of the simulation. It should be noted that the inference contains no approximations or arbitrary choices, and makes use of all the data. The key point of contention with observed data is whether the model is realistic to the system. If it is, then the inference faithfully converges on the underlying values regardless of what they are. If it is not, then this forms the point of departure for another round of modelling.

Figure 4.

Inference of the switching probabilities on simulated data given in Figure 1(c). The true value is shown as a dot in each plot. (Left) Joint marginal probability distributions for the exit probabilities for each state. (Right) Marginal probability distribution for each separate exit probability.

To illustrate the point, consider a model that is transitioning quickly between two states that are not very distinguishable. See for example the time series in Figure 5. It would be difficult to discern that there are two levels, let alone to estimate their switching rates. The model, and the assumption that there exists only two states, means that the inference will weigh each data point against the model and allow the extraction of the switching rates despite the noise.

Figure 5.

Time series with high relatively switching rates an difficult to discern states. Nevertheless the constraints imposed by the model allow the switching rate to be accurately estimated with enough data.

4. Inferring the state

A further inference we can make is to determine the underlying state behind a data point given all the data observed, but without knowledge of the transition probabilities $T_{\mathit{xy}}$ or the count rates $R_x$. Specifically, the task is to determine $P\left(s_k=x|\overrightarrow{c}I\right)$, where $x\in \left\{a,b,c,\dots \right\}$.

In order to do this inference we add $T,R$ and the rest of the states using marginalisation. An application of Bayes’ rule then yields,

$$\begin{array}{l}P\left(s_k=x|\overrightarrow{c}I\right)=\frac{1}{\mathcal{N}}\int \mathit{\text{dTdR}}\sum _{\overrightarrow{s}\ne s_x}P\left(\overrightarrow{c}|s_1\dots s_{k-1}\left(s_k=x\right)s_{k+1}\dots s_n\mathrm{\Omega }\right)\\ \times P\left(s_1\dots s_{k-1}\left(s_k=x\right)s_{k+1}\dots s_n|\mathrm{\Omega }\right)\end{array}$$E19

where $\mathrm{\Omega }=TRI$ as before, and the prior probability distributions where taken as constant with the normalisation $\mathcal{N}$ to be determined at the end. The summation is over all states $s_j$ where $j\ne k$ denoted in short as $\overrightarrow{s}\ne s_k$. Expanding the last term in temporal order and making use of the independencies in (3) and (4) leads to,

$$P\left(s_k=x|\overrightarrow{c}I\right)=\frac{1}{\mathcal{N}}\int \mathit{\text{dRdT}}\sum _{\overrightarrow{s}\ne s_k}\prod _{t=1}^{N}P\left(c_t{s}_t|s_{t-1}\mathrm{\Omega }\right)P\left(s_0|\mathrm{\Omega }\right).$$E20

The probability can again be efficiently computed using the matrices $Q_t$ (14) and vector ${\overrightarrow{D}}_0$ (15) but with a restricted choice for $Q_k$,

$$P\left(s_k=x|\overrightarrow{c}I\right)=\frac{1}{\mathcal{N}}\int \mathit{\text{dRdT}}{\overrightarrow{1}}^T{Q}_n\cdots Q_{k+1}{Q}_k^{\left(x\right)}{Q}_{k-1}\cdots Q_1{\overrightarrow{D}}_0.$$E21

where $Q_k^{\left(x\right)}$ only has the row corresponding to state $x$. So for the three-state example we would have

$$Q_k^{\left(b\right)}=\left[\begin{array}{ccc}0& 0& 0\\ P\left(c_k{s}_k=b|s_{k-1}=a\mathrm{\Omega }\right)& P\left(c_k{s}_k=b|s_{k-1}=b\mathrm{\Omega }\right)& P\left(c_k{s}_k=b|s_{k-1}=c\mathrm{\Omega }\right)\\ 0& 0& 0\end{array}\right].$$E22

5. Conclusions

We have derived efficient methods for obtaining the posterior distribution of transition rates given observation of accumulated signal for a hidden Markov model with multiple states. The methods do not require approximations to make the calculation tractable, and make use of all the observed data; the main fitting task is the selection and assumptions that go into the models. Moreover, the advantage of obtaining the posterior distribution is that rigorous error bounds can be placed on the inferred parameters for example by use of credible regions. We have also demonstrated how to use the data to determine the underlying state, and unlike the threshold technique, the result automatically carries error bounds in the form of a probability distribution. In general we expect the methods presented to apply to a wide variety of hidden Markov models with multiple states and other methods of indirect observation. The treatment here has looked at systems where there is a well-defined state time, which will apply to situations where the signal is sampled faster than the typical transition rates. These models are then capable of determining the model parameters for high or low switching probabilities. It is possible to generalise this approach to model systems which might switch at any time regardless of signal acquisition intervals [1].

Acknowledgments

This work was supported by the Australian Research Council (ARC) Centre of Excellence for Quantum Engineered Systems grant (CE170100009).

We recognise the intellectual and physical labour of this research was conducted on the traditional lands of the Wattamattagal clan of the Darug nation, and of the Pambalong clan of the Awabakal people.