Tuesday, March 31, 2009

Markov Process Sampling

One of the major strategies employed by single molecule techniques is to model particle traces with a Hidden Markov Model (HMM), assuming some number of states typically optimized by a Bayesian Information Criteria (BIC). These experiments are typically probing macroscopic fluctuations of molecular conformations by watching FRET (Forster Resonance Energy Transfer) trajectories. These trajectories are recorded, typically, with a fast CCD or perhaps a PMT (for non-wide-field studies). In the former case, integration times are on time scale of milliseconds.

I was curious about the nature of how these HMM react to undersampling. It's well known that undersampling a signal causes well known effects, like aliasing unless one takes care to sample below the Nyquist frequency. If the relevant underlying fluctuations of your system are faster than 500 microseconds, you're no longer above the Nyquist frequency for a CCD with a millisecond integration time. How does that affect the estimation of HMM parameters? It seems obvious to me that the transition matrix / emission matrix are going to be the only things affected, assuming you observe the system long enough that you see all possible states. But how exactly will the transition matrix be affected? Can it be compensated for?

It is pretty standard knowledge that undersampling causes the signal to be aliased to a lower frequency. It make intuitive sense to me that lower frequencies in HMMs are represented by lower transition probabilities, assuming the system has constant rates (which means that state duration can be modeled by an exponential distribution with expectated life time equal to the inverse of the transition probability). The reasoning here is that if the inverse of the state duration is the transition probability, then short state durations would result in high transition probabilities and long state durations will result in low transition probabilities.

To test this hypothesis, I decided to run a simulation in Matlab. I generated a 3-state sequence based off a known Markov Model with 1 million data points. I computed the upperbound of the 3 expected state durations and then went back and sampled the simulated sequence at a rate equal to 10 times the maximum expected state duration, to simulate under sampling. Sampling was periodic, with a small time variance for each point, and I took 300 data points for each "experiment". I then went back and estimated the transition probabilities for each of these "experiments" with a HMM algorithm. I averaged 50 of the resulting transition matrices. The results were surprising. Name, the exact opposite happened: the probabilities increased on undersampling. At certain degrees of undersampling (like 10x the maximum expected state duration), they seemed to converge to a multiple of the real transition probabilities related to the multiple of the inverse of the real transition probability. Other times, there didn't seem to be an obvious correlation.


I have yet to do much literature research on this particular problem, preferring to just jump in and see what I could come up with in an hour of coding. But it's obvious to me now that we're not going to be able to use standard signal analysis theory for sampling problems like this. For one, Markov sequences are non-periodic, so in retrospect, aliasing doesn't have much basis. I refuse (as of yet) to believe that one cannot predict the effects of undersampling (at least qualitatively), simply on the basis that no matter how I tweaked the parameters of my simulation the probabilities were always higher than the real ones, as long as I under sampled.


~~~~~~~~~~~~~Edit~~~~~~~~~~~~~~

On some additional thought, I think I can explain why the TM (transition matrix) is exhibiting an increase in probabilities. Since the sequence was generated from a first order Markov system, later time points are, by definition, not probabilistically linked to any point beyond the one previous to it, though they are weakly linked by conditional probabilities, which decrease rapidly in magnitude as one goes further. At some point, the coupling between observed events will be so weak that one may as well just be sampling a uniform distribution. I thought, originally, I was seeing a conservation of the relative ratios of the TM elements compared to the real TM. More experimentation proved this initial hope wrong. Here's a summary of those results:

1. Looking at the internally relative probabilistic ratios of the TMs, I am seeing fairly large fluctuation in how the observed and fit ratios compare to the ratios used to generate the data. It wasn't until I collected and averaged 1000 traces that I started to see reasonable convergence on the difference between observed and actual internal probability ratios in the TM. This has significant practical implications for experimental study of markov processes in a possibly undersampled regime.

2. I discovered a problem in my previous simulations such that the generating model's TM was ill defined (rows didn't sum to 1). I fixed this. The result was that the averaged, observed TM appears to be converging towards [0.3 0.3 0.3; 0.3 0.3 0.3; 0.3 0.3 0.3] (for a 3 state system), which is consistent with my earlier assertion that undersampling a Markov process might result in merely sampling a uniform distribution. I haven't tested how sampling rate affects this property yet, I've just made sure that I was under sampling the sequence such that I would sample on the time scale of ten times the off-diagonal transition characteristic time scale.

No comments:

Post a Comment