When a computer code is used to simulate a complex system, a fundamental task is to assess the uncertainty of the simulator. In the case of computationally expensive simulators, this is often accomplished via a surrogate statistical model, a statistical output emulator. An effective emulator is one that provides good approximations to the computer code output for wide ranges of input values. In addition, an emulator should be able to handle large dimensional simulation output for a relevant number of inputs; it should flexibly capture heterogeneities in the variability of the response surface; it should be fast to evaluate for arbitrary combinations of input parameters; and it should provide an accurate quantification of the emulation uncertainty.

In this work, we develop Bayesian adaptive spline methods for emulation of computer models that output functions. We introduce modifications to traditional Bayesian adaptive spline approaches that allow for fitting large amounts of data and allow for more efficient Markov chain Monte Carlo sampling. We develop a functional approach to sensitivity analysis that can be performed using this emulator. We present a sensitivity analysis of a computer model of the deformation of a protective plate used in pressure driven experiments. This example serves as an illustration of the ability of Bayesian adaptive spline emulators to fulfill all the necessities of computability, flexibility and reliable calculation on relevant measures of sensitivity.

We extend the methods to emulation of an atmospheric dispersion simulator that outputs a plume in space and time based on inputs detailing the characteristics of the release, some of which are categorical. We achieve accurate emulation using Bayesian adaptive splines to model weights on empirical orthogonal functions. We extend the adaptive spline methodology to allow for categorical inputs. We use this emulator as well as appropriately identifiable simulator discrepancy and observational error models to calibrate the simulator using a dataset from an experimental release of particles from the Diablo Canyon Nuclear Power Plant in Central California. Since the release was controlled, these characteristics are known, allowing us to compare our findings to the truth.

We further extend the methods to emulate a computer model that outputs misaligned functional data. We do this by modeling the aligned, or warped, data as well as the warping functions, using separate Bayesian adaptive spline models. We explore inference methods that treat these models jointly and separately, and establish methods to ensure that the warping functions are non-decreasing. These methods are applied to a high-energy-density physics model that outputs a curve representing energy as a function of time.

Computational efficiency is at the forefront of many cutting edge spatial modeling techniques. Non-stationarity, on the other hand, has often been an unintentional feature of the approximations used in these spatial models. Deriving from the well known multivariate linear regression model, we propose a non-stationary and non-isotropic spatial model. In order to remain relevant with today’s massive datasets challenges, we apply the concept of nearest-neighbors to our normal-inverse-Wishart framework. The model, called Nearest-Neighbor Gaussian Process with Random Covariance matrices (NN-RCM) is developed for both univariate and multivariate spatial settings and allow for specific characteristics such as duplicate observations and missing data. In a first dive into nearest-neighbor models such as Nearest-neighbor Gaussian Processes (NNGP), we initially look at divide-and-conquer implementations. Ultimately, we devise a technique to use NNGP models on blocks of data to then aggregate the posterior inference without loss of information. The models are illustrated in a case study of albedo assessments over CONUS from the Geostationary Operational Environmental Satellites (GOES) East and West. First, the divide-and-conquer algorithm is used to output multi-day successive predictions of albedo assessments in a sequential updating scheme. By applying the distributed model on each day of data, one can concatenate the posterior parameters of the blocks of interest to output evolving predictions. Finally, we apply the bivariate NN-RCM model using each satellite as a source of information. The objective is to merge the albedo assessments while also quantifying the discrepancy between the sources.

Large spatial datasets often exhibit fine scale features that only occur in sub-

domains of the space, coupled with large scale features at much larger ranges.

The most commonly used model used for spatial datasets is the Gaussian Process,

but evaluation of likelihood is computationally expensive. Additionally,

traditional Gaussian Processes models make very strong assumptions regarding

the symmetry of the Gaussian field. In particular they assume

stationarity, namely that covariance functions depend only on the displacement

vector between two points, not their locations. This assumption prevents stationary

Gaussian Processes from accounting for multi-scale features that only exist in parts

of the spatial domain. In this work, we develop multi-resolution kernel convolution methods that explicitly account for local multi-scale features through spatially varying resolution. These methods define an increasingly refined set of nested kernels, and induce sparsity on these grids.

We first introduce modifications to existing multi-resolution kernel convolution models

that result in spatially varying resolution through a sparsity inducing prior.

We cast spatially varying resolution as a model selection problem, and develop a

Shotgun Stochastic Search algorithm that considers an infinite number of resolutions,

and permits uncertainty quantification without resorting to MCMC.

We propose a LASSO like prior that achieves spatially varying resolution at its

maximum a posteriori, and develop a proximal gradient descent algorithm to find this

optimum considering an infinite number of resolutions. We then develop a Bayesian

model averaging approach to perform uncertainty quantification in this setting.

We implement these methods in an efficient and reproducible manner via the R package

MSSS. We discuss in detail the computational efficiency achieved by

leveraging parallel computation, compactly supported kernels, add one column

regression updates, and modern optimization methods in {\tt MSSS}. We demonstrate the local feature identification properties of spatially varying resolution and demonstrate the computational performance by considering a land surface temperature dataset from the Ozarks, and a large sea surface temperature dataset collected by a satellite off the coast of California

The focus of this thesis is to develop a general methodology to obtain high-resolution spatial-temporal forecasts of Sea Surface Temperature (SST) using an ensembles of general circulation model (GCM) output and historical records as the major driving force. As a case study, we consider Sea Surface Temperature (SST) in the North Pacific Ocean. We use two ensembles of different GCM simulation output, made available in the 4th Assessment Report of the Intergovernmental Panel on Climate Change: one corresponds to 20th century forcing conditions and the other corresponds to the A1B emissions scenario for the 21st century. Given a representation of the SST spatio-temporal fields based on a common set of empirical orthogonal functions (EOFs), we use a hierarchical Bayesian model for the EOF coefficients to estimate a baseline and a set of model discrepancies. These components are all time-varying. The model enables us to extract relevant temporal patterns of variability from both the observations and simulations and obtain common patterns from all eighteen series. This is used to obtain unified 21st century forecasts of relevant oceanic indexes as well as whole fields of forecast North Pacific SST. The unified forecast captures large longterm oceanic behavior, however the coarse resolution prevents us from capturing coastal behaviors. We use the unified forecast to model high resolution SST by establishing a link between large and small scale variability using statistical downscaling techniques. Using a combination of a discrete process convolution and a dynamic linear model, we obtain a smooth high-resolution forecast of SST fields off the coast of California. To model the high resolution data faster and efficiently, we developed and implement a parallel version of the forward filtering backwards sampling algorithm. We finish the work with remarks on the model results and address future avenues this work can take.

Model-based inferential methods for point processes have received less attention than the corresponding theory of point processes and is more scarcely developed than other areas of statistical inference.

Classical inferential methods for point processes include likelihood-based and nonparametric methods. Bayesian analysis provides simulation-based estimation of several statistics of interest for point processes. However, a challenge of Bayesian modeling, specifically for point processes, is selecting an appropriate parametric form for the intensity function. Bayesian nonparametric methods aim to avoid the narrow focus of parametric assumptions by imposing priors that can support the entire space of distributions and functions. It is naturally a more flexible and adaptable approach than those based on parametric models.

In this dissertation, we focus on developing methodology for some classes of temporal point processes modeling and inference in the context of Bayesian nonparametric methods, mainly with applications in environmental science.

Firstly, we are motivated to study seasonal marked point process by an application of

hurricanes occurrences. We develop nonparametric Bayesian methodology to study the dynamic

evolution of a seasonal marked point process intensity under the assumption that the point process is a non-homogeneous Poisson process. The dynamic model for time-varying intensities provides both the intra-seasonal and inter-seasonal variability of occurrences of events. Considering marks, we provide a full probabilistic model for the point process over the joint marks-points space which allows for different types of inferences, including full inference for dynamically evolving conditional mark densities given a time point, a particular time period, and even a subset of marks.

We apply this method to study the evolution of the intensity of the process of hurricane landfall occurrences, and the respective maximum wind speed and associated damages. We show several novel inferences which are explored for the first time in the analysis of hurricane occurrences.

Then we look beyond Poisson processes and propose a flexible approach to modeling and inference for homogeneous renewal processes.

This modeling method is based on a structured mixture of Erlang densities with common scale parameter for the renewal process inter-arrival density. The mixture weights are defined through an underlying distribution function modeled nonparametrically with a Dirichlet process prior. This model specification enables flexible shapes for the inter-arrival time density, including heavy tailed and multimodal densities. Moreover, the choice of the Dirichlet process centering distribution controls clustering or declustering patterns for the point process.

Finally we extend our model to accommodate point processes with time-varying inter-arrivals, which are referred to as modulated renewal processes in the literature. We introduce time dependence in the scale parameter of the Erlang mixture by replacing it with a latent stochastic process. A number of synthetic data sets and real data sets are used to illustrate the modeling approaches.

The main contribution of this thesis is to provide Bayesian nonparametric modeling and inference methods for some classes of point processes, which are more flexible than existing methods.

Moreover, the key complication for Bayesian inference is that the likelihood of a generic point process involves a normalizing constant which is, most of the times, analytically intractable. Discretization is very often used in existing methods to get likelihood approximations that facilitate computations, especially for models based on Gaussian process priors. Superior to these methods, our work uses the exact likelihood without approximation in all of our developed models.

Integro-Differential Equations (IDEs) are a novel way of dynamically modeling spatio-temporal data. IDEs are characterized by a kernel which controls the spatial and temporal associations. The ubiquitous choice for kernel has been Gaussian. We explore advantages of more flexible kernel choices. One-dimensional space is considered initially, replacing the Gaussian IDE kernel with more flexible parametric families of distributions. The kernels are chosen based on stochastic partial differential equation approximations which connect characteristics of the kernel with interpretable physical properties of the underlying process controlling the data. Next, Dirichlet process mixtures of normal distributions are used to model non-parametrically the IDE kernel. Computational issues arise using non-parametric kernels which are solved using Hermite polynomials and Hamiltonian Monte Carlo sampling. To develop flexible modeling in two-dimensional space, we propose bivariate stable distributions as IDE kernels. By using Bernstein polynomials as a prior for the measure defining the bivariate stable, a wide variety of shapes can be achieved. Bivariate stable kernels will be shown to outperform the Gaussian kernel by comparing K-step ahead predictions for Pacific sea surface temperature anomalies. Through study of properties for the proposed models, and empirical investigation with synthetic and real data, we demonstrate that the methodology has the potential to significantly improve the inference and forecasting capacity of IDE models based on Gaussian kernels.

This dissertation builds a modeling framework for non-Gaussian spatial processes, time series, and point processes, with a Bayesian inference paradigm that provides uncertainty quantification. Our methodological development emphasizes direct modeling of non-Gaussianity, in contrast with traditional approaches that consider data transformations or modeling through functionals of the data probability distribution. We achieve the goal by defining a joint distribution through factorization into a product of univariate conditional distributions according to a directed acyclic graph which implies conditional independence. We model each conditional distribution as a weighted combination of first-order conditionals, with weights that can be locally adaptive, for each one of a given number of parents which correspond to spatial nearest-neighbors or temporal lags. Such a formulation features specification of bivariate distributions that define the first-order conditionals for flexible, parsimonious modeling of multivariate non-Gaussian distributions. We obtain, in time, high-order Markov models with stationary marginals, and point process models for limited memory, dependent renewals, and duration clustering; and in space, nearest-neighbor mixture models for spatial processes. Regarding computation, representing the framework by directed acyclic graphs, with a mixture model formulation for the conditionals, gains efficiency and scalability relative to many non-Gaussian models. We develop Markov chain Monte Carlo algorithms for implementation of posterior inference and prediction, with data illustrations in biological, environmental, and social sciences.