Variance-reduced simulation of lattice Markov chains

The focus of this dissertation is on reducing the cost of Monte Carlo estimation for lattice-valued Markov chains. We achieve this goal by manipulating the random inputs to stochastic processes (Poisson random variables in the discrete-time setting and Poisson processes in continuous-time) such that they become negatively correlated with some of their cohort while their individual marginal distributions are completely unaltered. In this way, we preserve the convergence properties of the Law of Large Numbers, but mean estimates, say, constructed from these sample paths exhibit dramatically reduced variance. The work is comprised of three main parts. First, we introduce algorithms to reduce the simulation costs for discrete-time Markov chains. We describe how to modify the simulation of sample trajectories that introduces negative correlation while introducing no additional computational cost and that are compatible with existing codes. We support this algorithm with theoretical results, including guarantee that such mean estimators will be unbiased and consistent with respect to the discrete-time distribution. Further, we prove a recursive relation that characterizes the evolution of mutual negative covariance over time in the general case as well as prove a sufficient condition in the case of linear rate functions. Lastly, we present several numerical experiments that demonstrate multiple orders-of-magnitude reduction in mean-square error (MSE) for both linear and nonlinear reaction rate systems. In the next part, we show how insights gained from the discrete-time case can be used to inform a related approach in continuous-time. In these cases, we rely on a formulation of these lattice Markov chains called the random-time change representation. This allows us to translate the general problem of simulating anticorrelated trajectories of a given lattice Markov chain into the simpler problem of simulating anticorrelated pairs of unit-rate Poisson processes, which are the fundamental source of randomness that are input into random time-change representations. We systematically construct and analyze algorithms to produce negatively correlated, identically distributed Poisson processes. We prove closed form expressions for the MSE evolution of one of these systems, as well as present asymptotic performance lower bounds. We then show how to use these anticorrelated Poisson processes to simulate exact, identically distributed stochastic processes which are now significantly negatively correlated, and are thus suitable for variance-reduced Monte Carlo. Numerical experiments on both linear and nonlinear systems demonstrate order-of-magnitude cost reduction. We also introduce error vs cost comparisons with existing standard methods. Finally, we present extensions and refinements of the above algorithms. First is an approach to discrete- time simulation (specifically for tau-leaping systems) that leverages insights gained from the continuous-time approach in order to further strengthen the performance of the original algorithm in its weakest regime. This algorithm inherits several desirable properties from the antithetic discrete-time simulation case. In addition, we present numerical studies that show where this refinement outperforms the original algorithm. Finally, we present extensions of the anticorrelated simulation algorithms into both model predictive control and particle filtering

Applications of Markov Chain Monte Carlo methods to continuous gravitational wave data analysis

A new algorithm for the analysis of gravitational wave data from rapidly rotating neutron stars has been developed. The work is based on the Markov Chain Monte Carlo algorithm and features enhancements specifically targeted to this problem. The algorithm is tested on both synthetic data and hardware injections in the LIGO Hanford interferometer during its third science run ("S3''). By utilising the features of this probabilistic algorithm a search is performed for a rotating neutron star in the remnant of SN1987A within in frequency window of 4 Hz and a spindown window of 2E-10 Hz/s. A method for setting upper limits is described and used on this data in the absence of a detection setting an upper limit on strain of 7.3E-23. A further application of MCMC methods is made in the area of data analysis for the proposed LISA mission. An algorithm is developed to simultaneously estimate the number of sources and their parameters in a noisy data stream using reversible jump MCMC. An extension is made to estimate the position in the sky of a source and this is further improved by the implementation of a fast approximate calculation of the covariance matrix to enhance acceptance rates. This new algorithm is also tested upon synthetic data and the results are presented here. Conclusions are drawn from the results of this work, and comments are made on the development of MCMC algorithms within the field of gravitational wave data analysis, with a view to their increasing usage

Universal rank-order transform to extract signals from noisy data

We introduce an ordinate method for noisy data analysis, based solely on rank information and thus insensitive to outliers. The method is nonparametric and objective, and the required data processing is parsimonious. The main ingredients include a rank-order data matrix and its transform to a stable form, which provide linear trends in excellent agreement with least squares regression, despite the loss of magnitude information. A group symmetry orthogonal decomposition of the 2D rank-order transform for iid (white) noise is further ordered by principal component analysis. This two-step procedure provides a noise “etalon” used to characterize arbitrary stationary stochastic processes. The method readily distinguishes both the Ornstein-Uhlenbeck process and chaos generated by the logistic map from white noise. Ranking within randomness differs fundamentally from that in deterministic chaos and signals, thus forming the basis for signal detection. To further illustrate the breadth of applications, we apply this ordinate method to the canonical nonlinear parameter estimation problem of two-species radioactive decay, outperforming special-purpose least squares software. We demonstrate that the method excels when extracting trends in heavy-tailed noise and, unlike the Thiele-Sen estimator, is not limited to linear regression. A simple expression is given that yields a close approximation for signal extraction of an underlying, generally nonlinear signal

Doctor of Philosophy

dissertationFunctional magnetic resonance imaging (fMRI) measures the change of oxygen consumption level in the blood vessels of the human brain, hence indirectly detecting the neuronal activity. Resting-state fMRI (rs-fMRI) is used to identify the intrinsic functional patterns of the brain when there is no external stimulus. Accurate estimation of intrinsic activity is important for understanding the functional organization and dynamics of the brain, as well as differences in the functional networks of patients with mental disorders. This dissertation aims to robustly estimate the functional connectivities and networks of the human brain using rs-fMRI data of multiple subjects. We use Markov random field (MRF), an undirected graphical model to represent the statistical dependency among the functional network variables. Graphical models describe multivariate probability distributions that can be factorized and represented by a graph. By defining the nodes and the edges along with their weights according to our assumptions, we build soft constraints into the graph structure as prior information. We explore various approximate optimization methods including variational Bayesian, graph cuts, and Markov chain Monte Carlo sampling (MCMC). We develop the random field models to solve three related problems. In the first problem, the goal is to detect the pairwise connectivity between gray matter voxels in a rs-fMRI dataset of the single subject. We define a six-dimensional graph to represent our prior information that two voxels are more likely to be connected if their spatial neighbors are connected. The posterior mean of the connectivity variables are estimated by variational inference, also known as mean field theory in statistical physics. The proposed method proves to outperform the standard spatial smoothing and is able to detect finer patterns of brain activity. Our second work aims to identify multiple functional systems. We define a Potts model, a special case of MRF, on the network label variables, and define von Mises-Fisher distribution on the normalized fMRI signal. The inference is significantly more difficult than the binary classification in the previous problem. We use MCMC to draw samples from the posterior distribution of network labels. In the third application, we extend the graphical model to the multiple subject scenario. By building a graph including the network labels of both a group map and the subject label maps, we define a hierarchical model that has richer structure than the flat single-subject model, and captures the shared patterns as well as the variation among the subjects. All three solutions are data-driven Bayesian methods, which estimate model parameters from the data. The experiments show that by the regularization of MRF, the functional network maps we estimate are more accurate and more consistent across multiple sessions

Learning-accelerated algorithms for simulation and optimization

Simulation and optimization are fundamental building blocks for many computational methods in science and engineering. In this work, we explore the use of machine learning techniques to accelerate compute-intensive tasks in both simulation and optimization. Specifically, two algorithms are developed: (1) a variance reduction algorithm for Monte Carlo simulations of mean-field particle systems, and (2) a global optimization algorithm for noisy expensive functions. For the variance reduction algorithm, we develop an adaptive-control-variates technique for a class of simulations, where many particles interact via common mean fields. Due to the presence of a large number of particles and highly nonlinear dynamics, simulating these mean-field particle models is often time-consuming. Our algorithm treats the body of particles in the system as a source of training data, then uses machine learning to automatically build a model for the underlying particle dynamics, and finally constructs control variates with the learned model. We prove that the mean estimators from our algorithm are unbiased. More importantly, we show that, for a system with sufficiently many particles, our algorithm asymptotically produces more efficient estimators than naive Monte Carlo under certain regularity conditions. We applied our variance reduction algorithm to an aerosol particle simulation and found that the resulting simulation is about 7 times faster. The second algorithm is a parallel surrogate optimization algorithm, known as ProSRS, for noisy expensive black-box functions. Within this algorithm, we develop an efficient weighted-radial-basis regression procedure for constructing the surrogates. Furthermore, we introduce a novel tree-based technique, called the “zoom strategy”, to further improve optimization efficiency. We prove that if ProSRS is run for sufficiently long, with probability converging to one there will be at least one sample among all the evaluations that will be arbitrarily close to the global minimum. We compared ProSRS to several state-of-the-art Bayesian optimization algorithms on a suite of standard benchmark functions and two real machine-learning hyperparameter-tuning problems. We found that our algorithm not only achieves significantly faster optimization convergence, but is also orders of magnitude cheaper in computational cost. We also applied ProSRS to the problem of characterizing and validating a complex aerosol model against experimental measurements, where twelve simulation parameters must be optimized. This case illustrates the use of ProSRS for general global optimization problems

Mathematical Analysis and Modeling of Signaling Networks

Mathematical models are in focus of modern systems biology and increasingly important to understand and manipulate complex biological systems. At the same time, new and improved techniques in metabolomics and proteomics enhance the ability to measure cellular states and molecular concentrations. In consequence, this leads to important biological insights and novel potential drug targets. Model development in systems biology can be described as an iterative process of model refinement to match the observed properties. The resulting research cycle is based on a well-defined initial model and requires careful model revision in each step. {As an initial step, a stoichiometry-based mathematical model of the muscarinic acetylcholine receptor subtype 2 (M2 receptor)-induced signaling in Chinese hamster ovary (CHO) cells was derived. To validate the obtained initial model based on spatially accessible, not neces-sarily time-resolved data, the novel constrained flux sampling (CFS) is proposed in this work. The thus verified static model was then translated into a dynamical system based on ordinary differential equations (ODEs) by incorporating time-dependent experimental data. To learn from the errors of systems biological models, the dynamic elastic-net (DEN), a novel approach based on optimal control theory, is proposed in this thesis. Next, the Bayesian dy-namic elastic-net (BDEN), a systematic, fully algorithmic method based on the Markov chain Monte Carlo method was derived, which allows to detect hidden influences as well as missed reactions in ODE-based models. The BDEN allows for further validation of the developed M2 receptor-induced signaling pathway and thus provides evidence for the completeness of the obtained dynamical system. This thesis introduces the first comprehensive model of the M2 receptor-induced signaling in CHO cells. Furthermore, this work presents several novel algorithms to validate and correct static and dynamic models of biological systems in a semi-automatic manner. These novel algorithms are expected to simplify the development of further mathematical models in systems biology.Mathematische Modellierung und Analyse von Signalnetzwerken Mathematische Modelle stehen im Zentrum der modernen Systembiologie und werden immer wichtiger, um komplexe biologische Systeme verstehen und manipulieren zu können. Gleichzeitig erweitern neue und verbesserte Verfahren der Metabolomik und Proteomik die Möglichkeiten, Zellzustände und Molekülkonzentrationen zu bestimmen. Dies ermöglicht die Gewinnung neuer und wichtiger biologischer Erkenntnisse und die Identifizierung neuer potentieller Ansatzpunkte für medizinische Wirkstoffe. Die Modellentwicklung in der Systembiologie kann als ein iterativer Prozess der permanenten Modellverbesserung beschrieben werden, der das Ziel hat, die beobachteten Eigenschaften korrekt wiederzugeben. Der resultierende Modellierungskreislauf basiert auf einem klar bestimmten Anfangsmodell und erfordert das sorgfältige Anpassen des Modells in jedem einzelnen Modellierungsschritt. In einem ersten Schritt wurde ein auf stöchiometrischen Daten basierendes mathematisches Modell für die durch den muskarinischen Acetylcholinrezeptor des Subtyps 2 (M2-Rezeptor) induzierte Signalübertragung in CHO-Zellen aufgestellt. Zur Validierung des ursprünglichen Modells auf der Grundlage von räumlich erfassbaren, nicht notwendigerweise zeitaufgelösten Daten wird in dieser Arbeit das neu entwickelte Constrained Flux Sampling (CFS) vorgestellt. Das auf diese Weise verifizierte statische Modell wurde dann unter Einbeziehung zeitabhängiger experimenteller Messdaten in ein dynamisches Modell basierend auf gewöhnlichen Differentialgleichungen (DGL) umgewandelt. Um aus den mathematischen Unsicherheiten systembiologischer Modelle zu lernen, wird in dieser Arbeit das Dynamic Elastic-Net (DEN) eingeführt, ein neuer Ansatz basierend auf der Theorie der optimalen Steuerungen. Als nächster Schritt wurde das Bayesian Dynamic Elastic-Net (BDEN) entwickelt, eine systematische, vollständig algorithmische Methode basierend auf dem Markov-Chain-Monte-Carlo-Verfahren, die es erlaubt, sowohl verborgene Einflussfaktoren als auch übersehene Reaktionen in DGL-basierten Modellen aufzuspüren. Das BDEN ermöglicht die weitere Validierung des durch den M2-Rezeptor induzierten Signalwegs und liefert so den Beweis für die Vollständigkeit des modellierten dynamischen Systems. In dieser Arbeit wird das erste vollständige Modell für den durch den M2-Rezeptor induzierten Signalweg in CHO-Zellen eingeführt. Des Weiteren werden in dieser Arbeit verschiedene neue Algorithmen zur halbautomatischen Validierung und Korrektur statischer und dynamischer Modelle biologischer Systeme vorgestellt. Es wird erwartet, dass diese neuen Algorithmen die Entwicklung weiterer mathematischer Modelle in der Systembiologie stark vereinfachen

Modeling multivariate ultra-high-frequency financial data by Monte Carlo simulation methods

In questa tesi si propone una nuova classe di modelli probabilistici per dati multivariati ad altissima frequenza. Questi dati si incontrano oggigiorno in molti ambiti applicativi e in particolare in finanza quando si considerano contemporaneamente le transazioni di pi\uf9 di un\u2019azione. Le serie temporali di queste transazioni sono caratterizzate da tempi non equispaziati e non sincronizzati e per queste un naturale modello probabilistico di riferimento sono i processi puntuali marcati. In questo lavoro, per modellizzare questo tipo di dati abbiamo considerato una particolare sottoclasse di questi processi, in particolare la classe dei processi di Poisson doppio stocastici con marchi. Nello specifico del caso multivariato, si \ue8 assunto, per ogni azione, che i tempi di arrivo delle transazioni fossero descrivibili da un processo di Poisson doppio stocastico e che le relative intensit\ue0 latenti fossero funzione di alcune componenti dinamiche specifiche e di una componente dinamica comune, tutte di forma \u201cshot noise\u201d. Quest\u2019ultima componente dovrebbe essere responsabile del comportamento osservato sul mercato di alcuni panieri di azioni. Il problema principale posto da questa classe di modelli \ue8 il filtraggio delle intensit\ue0 latenti non osservabili sulla base delle transazioni osservate. Nella tesi si \ue8 proposto di affrontare questo problema di filtraggio non lineare ideando ed implementando una procedura stocastica basata sull\u2019algoritmo \u201creversibile jump Markov chain Monte Carlo\u201d. Per mezzo di questo algoritmo, si \ue8 riusciti a ricostruire a posteriori, non solo le intensit\ue0 latenti, ma anche le loro componenti, in particolare quella comune. Da un punto di vista empirico, sulla base di innumerevoli confronti tra le propriet\ue0 statistiche, relative principalmente alle correlazioni e alle cross-correlazioni tra coppie di azioni, di dati reali provenienti dalla Borsa di Milano e dati simulati, ottenuti sulla base di diverse ipotesi per i tempi di arrivo delle transazioni e per i logreturns, il modello proposto \ue8 risultato essere il pi\uf9 plausibile fornendo quindi un\u2019evidenza empirica per l\u2019esistenza di una componente comune sottostante i tempi di arrivo delle transazioni di panieri di azioni.In this thesis, we propose a modeling framework for multivariate ultra-high-frequency financial data. The proposed models belong to the class of the doubly stochastic Poisson processes with marks which are characterized by the number of events in any time interval to be conditionally Poisson distributed, given another positive stochastic process called intensity. The key assumption of these models is that the intensities are specified through a latent common dynamic factor that jointly drives their common behavior. Assuming the intensities are unobservable, we propose a signal extraction (filtering) method based on the reversible jump Markov chain Monte Carlo algorithm. Our proposed filtering method allows to filter not only the intensities but also their specific and common components. From an empirical stand point, on the basis of a comparison of real data with Monte Carlo simulated data, obtained under different assumptions for ticks (times and logreturns), based mainly on the behavior of the correlation between pairs of assets as a function of the sampling period (Epps effect), we found evidence for the existence of a single latent common factor responsible for the behavior observed in a set of assets from the Borsa di Milano