Temporal data modeling plays a vital role in various research including finance, environmental science and neuroscience. Understanding and interpreting the evolutionary system behind temporal data is of interest. This work mainly emphasizes efficient statistical models on temporal data via stochastic processes. In particular, we focus on statistical modeling via two flexible random processes: Markov processes and Gaussian processes.

The first stage of the research involves a novel hidden Markov model based on Markov jump processes for cervical cancer screening test data. This model is able to model the heterogeneity of both individual and time. We provide an efficient and scalable expectation maximization based inference approach. To the best of our knowledge, our model is the first statistical model that is able to scale to a population-level dataset.

Next, we consider an alternative stochastic process that is widely applied in temporal data modeling, the Gaussian process. Motivated by the kernel perspective of inducing-point based sparse Gaussian processes, we propose a general regularization framework of sparse Gaussian processes and extend it into latent variable models. We specifically consider variational inference under our regularization framework in various settings. We theoretically demonstrate that the variational inference under our regularization can be treated by maximizing a log likelihood lower bound on a corresponding empirical Bayesian model. Our framework is illustrated in various settings throughout both synthetic and real datasets.

Building on our proposed regularization framework, we develop a hierarchical sparse latent Gaussian process model specifically for categorical data and then we extend our model to temporal data via dynamical priors. In particular, we propose efficient variational inference to make it applicable to large datasets. Moreover, our model provides a visualization way to summarize the dynamics of categorical data into a low-dimension manifold.

The fourth project is motivated by electronic health record data and spatially varying coefficient linear coregionalization model. We propose a novel nonstationary multivariate Gaussian process model that allows it to model time dependent smoothness, scale and correlation across different dimensions. One special case of our model is emphasized due to its computational efficiency. It allows an efficient inference via Kronecker algebra. Moreover, we provide both Hamiltonian Monte Carlo inference as a fully Bayesian inference and Maximum a Posteriori based inference as an approximate Bayesian inference. Our posterior inference provides a promising way to understand the relation between cross-correlation of clinical variables and health status, which contributes to early disease detection.

Expensive black box systems arise in many engineering applications but can be difficult to optimize because their output functions may be complex, multi-modal, and difficult to understand. The task becomes even more challenging when the optimization is subject to constraints and no derivative information is available. In this dissertation, we combine response surface modeling, sequential Monte Carlo, and filter methods in order to solve problems of this nature. Furthermore, we propose a new model for correlated outputs of mixed type. Our modeling framework extends Gaussian process methodology for modeling of continuous multivariate spatial outputs by adding a latent process structure that allows for joint modeling of a variety of types of correlated outputs. In addition, we implement fully Bayesian inference using particle learning, which allows us to conduct fast sequential inference. By employing a filter algorithm for solving constrained optimization problems, we also establish two novel probabilistic metrics for guiding the filter. We extend these ideas to a multidimensional filter that outperforms the traditional filter method. Overall, this hybridization of statistical modeling and nonlinear programming efficiently utilizes both global and local search in order to converge to a global solution to the constrained optimization problem. To demonstrate the effectiveness of the proposed methods, we perform numerical tests on both synthetic examples and real-world hydrology computer experiment optimization problems.

The general problem addressed in this work concerns the analysis of a function of multiple real variables in which the output is itself a function of a real variable as well as categorical information. We are interested in the prediction of the output and the analysis of the boundaries separating regions during classification. Efficient and accurate prediction of the output curves is important for safety analysis and validation of a complex system, like an adaptive flight control system. An understanding of the boundaries between regions, where the aircraft can maintain stable flight or will break up is essential for safety certification. As a motivating application we are using the NASA IFCS (Intelligent Flight Control System) flight simulator.

For output prediction, we developed a new statistical method for emulation of computer models with multiple outputs. An emulator is a computationally efficient statistical model that is used to approximate a computationally expensive simulation by treating the simulator as a black box and learning a mapping from inputs to outputs. Our approach for emulating the curves is to represent them in an orthogonal basis (e.g., PCA, Fourier, or Wavelet), which captures curve characteristics, and then to predict the coefficients. We allow for the possibility of output curves whose length varies with input, which may occur when a simulator fails to run to completion for some inputs and the failures occur at different output time steps. To the best of our knowledge, the variable-length output problem has not yet been addressed. We have developed a hierarchical model, which first uses classification into a few groups and then fitting distinct models for different classes of output curves.

For the analysis of boundaries, we developed a new sequential approach based upon design of computer experiments. A dictionary of suitable linear or non-linear parameterized boundary shapes, which capture underlying physical and design knowledge can be provided by the domain expert. We incorporate this knowledge into our modeling and determine the most likely shapes and its parameters. Since each iteration requires a costly run of the system simulator, we developed a candidate selection mechanism, which is specifically tailored toward boundary detection and which can take priors into account in order to reduce the number of required simulation runs. We present results of experiments with artificial and IFCS simulation data sets.

Fisheries stock assessments often assume a two-parameter functional form (e.g. Beverton-Holt, Ricker, or Logistic) for the expected recruitment produced by a given level of spawning output. Mangel et. al. (2013) and others have shown that biological reference points (RP) such as Fmsy/M and Bmsy/B0 are largely determined by a single parameter (steepness) when using two-parameter relationships. These functions introduce strong correlations between RPs that are pre-determined by the functional form, rather than a biological characteristic of the stock. Mangel et. al. note that use of a three-parameter stock-recruitment relationship allows for independent estimation of these RPs. Built around these ideas, a novel simulation framework is developed to investigate the nature of biases in RP estimates that results from fitting a two-parameter functional form when the true relationship is more complicated. First methods for generating space-filling simulation designs in the RP space of three-parameter models are developed. By simulating misspecified RP inference under common two-parameter models over these simulation designs a Gaussian Process metamodel of two-parameter RP inference is developed to control for a spectrum of common ways that two-parameter models are misspecified. This analysis demonstrates the useful limits of commonly used population dynamics models, informs the utility of reducing RP bias, and suggests mechanisms for understanding how, and when, the most common two-parameter models fail to estimate RPs. The studied models vary in complexity from the Schaefer model to delay differential models including dynamics of individual growth and lagged maturity. Additionally, the methods presented can easily be extended to further include age-based frameworks.

In nanopore experiments, DNA replication facilitated by $\phi 29$ DNA polymerase (DNAP) can be observed at the single molecule level. The biochemical state of the DNA-DNAP complex was studied by setting the complex atop a $\alpha$-hemolysin nanopore and applying an electric voltage. The movement of the DNA strand relative to the nanopore was observed on a single base pair level by the ionic current blockade. The time trace of the recorded ionic current amplitude from these experiments was used to study the biochemical states. Given that the recorded amplitude of the ionic current was an indirect measurement of the true amplitude level, which in turn was an indirect measurement of the true biochemical state, the experiments were modeled as a Hidden Markov Chain (HMC). When the DNA position of two biochemical states relative to the nanopore is the same, the states yield the same current amplitude level. To extract the dynamic transition rates between biochemical states that are not distinguishable in amplitude level, two methodologies were applied to study the HMC. The first was a fully Bayesian model, for which Markov chain Monte-Carlo (MCMC) simulations were used to infer the reaction rates in a system of three biochemical states with two observed amplitude levels. The second model adopted concepts of Viterbi training or the segmental k-means algorithm to find point estimates of the transition rates. Given the low transition probabilities, the properties of the second model led to a substantial bias in inference. The bias was addressed by first using a meta-model to describe the relation between the generating transition rates and the biased inference. Then the inverse problem of the meta-model was solved to reduce the bias in the inference. The meta-model was a fully Bayesian Gaussian process model, built by creating a series of computer-generated datasets around the given dataset. Improved inference was obtained by drawing posterior samples from the inverse of the meta-model. Since both the MCMC simulations and bias reduction techniques resulted in simulated posteriors, the two methods were used to confirm each other. In the comparison of these two methods, the approach of the second model plus bias reduction has the advantage of achieving similar inference accuracy at a much lower computational cost.