Hidden state models for time series

Amongst all the objectives in the study of time series, uncovering the dynamic law of its generation is probably the most important. When the underlying dynamics are not available, time series modelling consists of developing a model which best explains a sequence of observations. In this thesis, we consider hidden space models for analysing and describing time series. We first provide an introduction to the principal concepts of hidden state models and draw an analogy between hidden Markov models and state space models. Central ideas such as hidden state inference or parameter estimation are reviewed in detail. A key part of multivariate time series analysis is identifying the delay between different variables. We present a novel approach for time delay estimating in a non-stationary environment. The technique makes use of hidden Markov models and we demonstrate its application for estimating a crucial parameter in the oil industry. We then focus on hybrid models that we call dynamical local models. These models combine and generalise hidden Markov models and state space models. Probabilistic inference is unfortunately computationally intractable and we show how to make use of variational techniques for approximating the posterior distribution over the hidden state variables. Experimental simulations on synthetic and real-world data demonstrate the application of dynamical local models for segmenting a time series into regimes and providing predictive distributions.

A pragmatic Bayesian approach to wind field retrieval

The ERS-1 Satellite was launched in July 1991 by the European Space Agency into a polar orbit at about 800 km, carrying a C-band scatterometer. A scatterometer measures the amount of backscatter microwave radiation reflected by small ripples on the ocean surface induced by sea-surface winds, and so provides instantaneous snap-shots of wind flow over large areas of the ocean surface, known as wind fields. Inherent in the physics of the observation process is an ambiguity in wind direction; the scatterometer cannot distinguish if the wind is blowing toward or away from the sensor device. This ambiguity implies that there is a one-to-many mapping between scatterometer data and wind direction. Current operational methods for wind field retrieval are based on the retrieval of wind vectors from satellite scatterometer data, followed by a disambiguation and filtering process that is reliant on numerical weather prediction models. The wind vectors are retrieved by the local inversion of a forward model, mapping scatterometer observations to wind vectors, and minimising a cost function in scatterometer measurement space. This thesis applies a pragmatic Bayesian solution to the problem. The likelihood is a combination of conditional probability distributions for the local wind vectors given the scatterometer data. The prior distribution is a vector Gaussian process that provides the geophysical consistency for the wind field. The wind vectors are retrieved directly from the scatterometer data by using mixture density networks, a principled method to model multi-modal conditional probability density functions. The complexity of the mapping and the structure of the conditional probability density function are investigated. A hybrid mixture density network, that incorporates the knowledge that the conditional probability distribution of the observation process is predominantly bi-modal, is developed. The optimal model, which generalises across a swathe of scatterometer readings, is better on key performance measures than the current operational model. Wind field retrieval is approached from three perspectives. The first is a non-autonomous method that confirms the validity of the model by retrieving the correct wind field 99% of the time from a test set of 575 wind fields. The second technique takes the maximum a posteriori probability wind field retrieved from the posterior distribution as the prediction. For the third technique, Markov Chain Monte Carlo (MCMC) techniques were employed to estimate the mass associated with significant modes of the posterior distribution, and make predictions based on the mode with the greatest mass associated with it. General methods for sampling from multi-modal distributions were benchmarked against a specific MCMC transition kernel designed for this problem. It was shown that the general methods were unsuitable for this application due to computational expense. On a test set of 100 wind fields the MAP estimate correctly retrieved 72 wind fields, whilst the sampling method correctly retrieved 73 wind fields.

Bayesian data assimilation

This thesis addresses data assimilation, which typically refers to the estimation of the state of a physical system given a model and observations, and its application to short-term precipitation forecasting. A general introduction to data assimilation is given, both from a deterministic and' stochastic point of view. Data assimilation algorithms are reviewed, in the static case (when no dynamics are involved), then in the dynamic case. A double experiment on two non-linear models, the Lorenz 63 and the Lorenz 96 models, is run and the comparative performance of the methods is discussed in terms of quality of the assimilation, robustness "in the non-linear regime and computational time. Following the general review and analysis, data assimilation is discussed in the particular context of very short-term rainfall forecasting (nowcasting) using radar images. An extended Bayesian precipitation nowcasting model is introduced. The model is stochastic in nature and relies on the spatial decomposition of the rainfall field into rain "cells". Radar observations are assimilated using a Variational Bayesian method in which the true posterior distribution of the parameters is approximated by a more tractable distribution. The motion of the cells is captured by a 20 Gaussian process. The model is tested on two precipitation events, the first dominated by convective showers, the second by precipitation fronts. Several deterministic and probabilistic validation methods are applied and the model is shown to retain reasonable prediction skill at up to 3 hours lead time. Extensions to the model are discussed.

A Note on Bayesian Estimation for the Negative-Binomial Model

2000 Mathematics Subject Classification: 62F15.

Variational Markov chain Monte Carlo for Bayesian smoothing of non-linear diffusions

In this paper we develop set of novel Markov Chain Monte Carlo algorithms for Bayesian smoothing of partially observed non-linear diffusion processes. The sampling algorithms developed herein use a deterministic approximation to the posterior distribution over paths as the proposal distribution for a mixture of an independence and a random walk sampler. The approximating distribution is sampled by simulating an optimized time-dependent linear diffusion process derived from the recently developed variational Gaussian process approximation method. The novel diffusion bridge proposal derived from the variational approximation allows the use of a flexible blocking strategy that further improves mixing, and thus the efficiency, of the sampling algorithms. The algorithms are tested on two diffusion processes: one with double-well potential drift and another with SINE drift. The new algorithm's accuracy and efficiency is compared with state-of-the-art hybrid Monte Carlo based path sampling. It is shown that in practical, finite sample applications the algorithm is accurate except in the presence of large observation errors and low to a multi-modal structure in the posterior distribution over paths. More importantly, the variational approximation assisted sampling algorithm outperforms hybrid Monte Carlo in terms of computational efficiency, except when the diffusion process is densely observed with small errors in which case both algorithms are equally efficient. © 2011 Springer-Verlag.

Bayesian Inference in Large-scale Problems

Many modern applications fall into the category of "large-scale" statistical problems, in which both the number of observations n and the number of features or parameters p may be large. Many existing methods focus on point estimation, despite the continued relevance of uncertainty quantification in the sciences, where the number of parameters to estimate often exceeds the sample size, despite huge increases in the value of n typically seen in many fields. Thus, the tendency in some areas of industry to dispense with traditional statistical analysis on the basis that "n=all" is of little relevance outside of certain narrow applications. The main result of the Big Data revolution in most fields has instead been to make computation much harder without reducing the importance of uncertainty quantification. Bayesian methods excel at uncertainty quantification, but often scale poorly relative to alternatives. This conflict between the statistical advantages of Bayesian procedures and their substantial computational disadvantages is perhaps the greatest challenge facing modern Bayesian statistics, and is the primary motivation for the work presented here.

Two general strategies for scaling Bayesian inference are considered. The first is the development of methods that lend themselves to faster computation, and the second is design and characterization of computational algorithms that scale better in n or p. In the first instance, the focus is on joint inference outside of the standard problem of multivariate continuous data that has been a major focus of previous theoretical work in this area. In the second area, we pursue strategies for improving the speed of Markov chain Monte Carlo algorithms, and characterizing their performance in large-scale settings. Throughout, the focus is on rigorous theoretical evaluation combined with empirical demonstrations of performance and concordance with the theory.

One topic we consider is modeling the joint distribution of multivariate categorical data, often summarized in a contingency table. Contingency table analysis routinely relies on log-linear models, with latent structure analysis providing a common alternative. Latent structure models lead to a reduced rank tensor factorization of the probability mass function for multivariate categorical data, while log-linear models achieve dimensionality reduction through sparsity. Little is known about the relationship between these notions of dimensionality reduction in the two paradigms. In Chapter 2, we derive several results relating the support of a log-linear model to nonnegative ranks of the associated probability tensor. Motivated by these findings, we propose a new collapsed Tucker class of tensor decompositions, which bridge existing PARAFAC and Tucker decompositions, providing a more flexible framework for parsimoniously characterizing multivariate categorical data. Taking a Bayesian approach to inference, we illustrate empirical advantages of the new decompositions.

Latent class models for the joint distribution of multivariate categorical, such as the PARAFAC decomposition, data play an important role in the analysis of population structure. In this context, the number of latent classes is interpreted as the number of genetically distinct subpopulations of an organism, an important factor in the analysis of evolutionary processes and conservation status. Existing methods focus on point estimates of the number of subpopulations, and lack robust uncertainty quantification. Moreover, whether the number of latent classes in these models is even an identified parameter is an open question. In Chapter 3, we show that when the model is properly specified, the correct number of subpopulations can be recovered almost surely. We then propose an alternative method for estimating the number of latent subpopulations that provides good quantification of uncertainty, and provide a simple procedure for verifying that the proposed method is consistent for the number of subpopulations. The performance of the model in estimating the number of subpopulations and other common population structure inference problems is assessed in simulations and a real data application.

In contingency table analysis, sparse data is frequently encountered for even modest numbers of variables, resulting in non-existence of maximum likelihood estimates. A common solution is to obtain regularized estimates of the parameters of a log-linear model. Bayesian methods provide a coherent approach to regularization, but are often computationally intensive. Conjugate priors ease computational demands, but the conjugate Diaconis--Ylvisaker priors for the parameters of log-linear models do not give rise to closed form credible regions, complicating posterior inference. In Chapter 4 we derive the optimal Gaussian approximation to the posterior for log-linear models with Diaconis--Ylvisaker priors, and provide convergence rate and finite-sample bounds for the Kullback-Leibler divergence between the exact posterior and the optimal Gaussian approximation. We demonstrate empirically in simulations and a real data application that the approximation is highly accurate, even in relatively small samples. The proposed approximation provides a computationally scalable and principled approach to regularized estimation and approximate Bayesian inference for log-linear models.

Another challenging and somewhat non-standard joint modeling problem is inference on tail dependence in stochastic processes. In applications where extreme dependence is of interest, data are almost always time-indexed. Existing methods for inference and modeling in this setting often cluster extreme events or choose window sizes with the goal of preserving temporal information. In Chapter 5, we propose an alternative paradigm for inference on tail dependence in stochastic processes with arbitrary temporal dependence structure in the extremes, based on the idea that the information on strength of tail dependence and the temporal structure in this dependence are both encoded in waiting times between exceedances of high thresholds. We construct a class of time-indexed stochastic processes with tail dependence obtained by endowing the support points in de Haan's spectral representation of max-stable processes with velocities and lifetimes. We extend Smith's model to these max-stable velocity processes and obtain the distribution of waiting times between extreme events at multiple locations. Motivated by this result, a new definition of tail dependence is proposed that is a function of the distribution of waiting times between threshold exceedances, and an inferential framework is constructed for estimating the strength of extremal dependence and quantifying uncertainty in this paradigm. The method is applied to climatological, financial, and electrophysiology data.

The remainder of this thesis focuses on posterior computation by Markov chain Monte Carlo. The Markov Chain Monte Carlo method is the dominant paradigm for posterior computation in Bayesian analysis. It has long been common to control computation time by making approximations to the Markov transition kernel. Comparatively little attention has been paid to convergence and estimation error in these approximating Markov Chains. In Chapter 6, we propose a framework for assessing when to use approximations in MCMC algorithms, and how much error in the transition kernel should be tolerated to obtain optimal estimation performance with respect to a specified loss function and computational budget. The results require only ergodicity of the exact kernel and control of the kernel approximation accuracy. The theoretical framework is applied to approximations based on random subsets of data, low-rank approximations of Gaussian processes, and a novel approximating Markov chain for discrete mixture models.

Data augmentation Gibbs samplers are arguably the most popular class of algorithm for approximately sampling from the posterior distribution for the parameters of generalized linear models. The truncated Normal and Polya-Gamma data augmentation samplers are standard examples for probit and logit links, respectively. Motivated by an important problem in quantitative advertising, in Chapter 7 we consider the application of these algorithms to modeling rare events. We show that when the sample size is large but the observed number of successes is small, these data augmentation samplers mix very slowly, with a spectral gap that converges to zero at a rate at least proportional to the reciprocal of the square root of the sample size up to a log factor. In simulation studies, moderate sample sizes result in high autocorrelations and small effective sample sizes. Similar empirical results are observed for related data augmentation samplers for multinomial logit and probit models. When applied to a real quantitative advertising dataset, the data augmentation samplers mix very poorly. Conversely, Hamiltonian Monte Carlo and a type of independence chain Metropolis algorithm show good mixing on the same dataset.

Mixtures of g-priors in Generalized Linear Models

Mixtures of Zellner's g-priors have been studied extensively in linear models and have been shown to have numerous desirable properties for Bayesian variable selection and model averaging. Several extensions of g-priors to Generalized Linear Models (GLMs) have been proposed in the literature; however, the choice of prior distribution of g and resulting properties for inference have received considerably less attention. In this paper, we extend mixtures of g-priors to GLMs by assigning the truncated Compound Confluent Hypergeometric (tCCH) distribution to 1/(1+g) and illustrate how this prior distribution encompasses several special cases of mixtures of g-priors in the literature, such as the Hyper-g, truncated Gamma, Beta-prime, and the Robust prior. Under an integrated Laplace approximation to the likelihood, the posterior distribution of 1/(1+g) is in turn a tCCH distribution, and approximate marginal likelihoods are thus available analytically. We discuss the local geometric properties of the g-prior in GLMs and show that specific choices of the hyper-parameters satisfy the various desiderata for model selection proposed by Bayarri et al, such as asymptotic model selection consistency, information consistency, intrinsic consistency, and measurement invariance. We also illustrate inference using these priors and contrast them to others in the literature via simulation and real examples.

Haplotype distribution of five nuclear genes based on network genealogies and Bayesian inference indicates that Trypanosoma cruzi hybrid strains are polyphyletic

Chagas disease is still a major public health problem in Latin America. Its causative agent, Trypanosoma cruzi, can be typed into three major groups, T. cruzi I, T. cruzi II and hybrids. These groups each have specific genetic characteristics and epidemiological distributions. Several highly virulent strains are found in the hybrid group; their origin is still a matter of debate. The null hypothesis is that the hybrids are of polyphyletic origin, evolving independently from various hybridization events. The alternative hypothesis is that all extant hybrid strains originated from a single hybridization event. We sequenced both alleles of genes encoding EF-1 alpha, actin and SSU rDNA of 26 T. cruzi strains and DHFR-TS and TR of 12 strains. This information was used for network genealogy analysis and Bayesian phylogenies. We found T. cruzi I and T. cruzi II to be monophyletic and that all hybrids had different combinations of T. cruzi I and T. cruzi II haplotypes plus hybrid-specific haplotypes. Bootstrap values (networks) and posterior probabilities (Bayesian phylogenies) of clades supporting the monophyly of hybrids were far below the 95% confidence interval, indicating that the hybrid group is polyphyletic. We hypothesize that T. cruzi I and T. cruzi II are two different species and that the hybrids are extant representatives of independent events of genome hybridization, which sporadically have sufficient fitness to impact on the epidemiology of Chagas disease.

Distribution of an asymmetrical copepod, Hatschekia plectropomi, on the gills of Plectropomus leopardus

Hatschekia plectropomi, an ectoparasitic copepod found on the gills, infected Plectropomus leopardus from Heron Island Reef with 100% prevalence (n = 32) and a mean +/- S.E. infection intensity of 131.9 +/- 22.1. The distribution of 4222 adult female parasites across 32 individual host fish was investigated at several organizational levels ranging from the level of holobranch pairs to that of individual filaments. Parasites demonstrated a site preference for the two central holobranchs (2 and 3). Along the lengths of hemibranchs, filaments near the dorsal and ventral ends and those in the proximity of the bend region were rarely occupied. The probability of coming into contact with a suitable attachment site and the ability to withstand ventilation forces at that site were proposed as the major factors affecting distribution. Two H. plectropomi morphotypes were identified based on the direction of body curvature. Regardless of morphotype, 99.9% of individuals were attached such that the convex side of the body was oriented towards the oncoming ventilating water currents. Further, 93.3% of individuals attached to the posterior faces of filaments, leading to a predictable pattern of attachment for this species. It is suggested that the direction of body curvature develops in response to the direction of the ventilating water currents. (c) 2006 The Fisheries Society of the British Isles.

Distribution of lymph nodes in the mesorectum: how deep is TME necessary?

Background Standardization of total mesorectal excision (TME) had a great impact on decreasing local recurrence rates for the treatment of rectal cancer. However, exact numbers and distribution of lymph nodes (LN) along the mesorectum remains controversial with some studies suggesting that few LNs are present in the distal third of the mesorectum. Methods Eighteen fresh cadavers without a history of rectal cancer were studied. The rectum was removed by TME and then was divided into right lateral, posterior and left lateral sides, which were further subdivided into 3 levels (upper, middle and lower). A pathologist determined the number and sizes of the LNs in each of the nine areas, b linded to their anatomical origin. Results Overall, the mesorectum had a mean of 5.7 LNs (SD=3.7) and on average each LN had a maximum diameter of 3.0 mm (SD=2.7). There was no association between the mean number or size of LNs with gender, BMI, or age. There was a significantly higher prevalence of LNs in the posterior location (2.8 per mesorectum) than in the two lateral locations (0.8 and 1.2 per mesorectum; p=0.02). The distribution of LNs in the three levels of the rectum was not significant. Conclusions The distribution of LNs reinforces the fact that TME should always include the distal third of the mesorectum. Care must be taken to not violate the posterior aspect of the mesorectum.

Distribution kinetics of solutes in the isolated in-situ perfused rat head using the multiple indicator dilution technique and a physiological two-barrier model

The purpose of this study was to determine the pharmacokinetics of [C-14]diclofenac, [C-14]salicylate and [H-3]clonidine using a single pass rat head perfusion preparation. The head was perfused with 3-[N-morpholino] propane-sulfonic acid-buffered Ringer's solution. Tc-99m-red blood cells and a drug were injected in a bolus into the internal carotid artery and collected from the posterior facial vein over 28 min. A two-barrier stochastic organ model was used to estimate the statistical moments of the solutes. Plasma, interstitial and cellular distribution volumes for the solutes ranged from 1.0 mL (diclofenac) to 1.6 mL (salicylate), 2.0 mL (diclofenac) to 4.2 mL (water) and 3.9 mL (salicylate) to 20.9 mL (diclofenac), respectively. A comparison of these volumes to water indicated some exclusion of the drugs from the interstitial space and salicylate from the cellular space. Permeability-surface area (PS) products calculated from plasma to interstitial fluid permeation clearances (CLPI) (range 0.02-0.40 mL s(-1)) and fractions of solute unbound in the perfusate were in the order: diclofenac>salicylate >clonidine>sucrose (from 41.8 to 0.10 mL s(-1)). The slow efflux of diclofenac, compared with clonidine and salicylate, may be related to its low average unbound fraction in the cells. This work accounts for the tail of disposition curves in describing pharmacokinetics in the head.

The distribution and morphological characteristics of cholinergic cells in the brain of monotremes as revealed by ChAT immunohistochemistry

The present study employs choline acetyltransferase (ChAT) immunohistochemistry to identify the cholinergic neuronal population in the central nervous system of the monotremes. Two of the three extant species of monotreme were studied: the platypus (Omithorhynchus anatinus) and the short-beaked echidna (Tachyglossus aculeatus). The distribution of cholinergic cells in the brain of these two species was virtually identical. Distinct groups of cholinergic cells were observed in the striatum, basal forebrain, habenula, pontomesencephalon, cranial nerve motor nuclei, and spinal cord. In contrast to other tetrapods studied with this technique, we failed to find evidence for cholinergic cells in the hypothalamus, the parabigeminal nucleus (or nucleus isthmus), or the cerebral cortex. The lack of hypothalamic cholinergic neurons creates a hiatus in the continuous antero-posterior aggregation of cholinergic neurons seen in other tetrapods. This hiatus might be functionally related to the phenomenology of monotreme sleep and to the ontogeny of sleep in mammals, as juvenile placental mammals exhibit a similar combination of sleep elements to that found in adult monotremes. Copyright (C) 2002 S. Karger AG, Basel.

Surgical revascularization of posterior coronary arteries without cardiopulomonary bypass

OBJECTIVE: To assess the results observed during the early postoperative period in patients who had the posterior coronary arteries revascularized without cardiopulmonary bypass (CPB), in regard to the following parameters: age, sex,bypass grafts types, morbidity and mortality. METHODS: From January 1995 to June 1998, 673 patients underwent myocardial revascularization (MR). Of this total, 607 (90.20%) MR procedures were performed without CPB. The posterior coronary arteries (PCA) were revascularized in 298 (44.27%) patients, 280 (93.95%) without CPB. The age of the patients ranged from 37 to 88 years (mean, 61 years). The male gender predominated, with 198 men (70.7%). The revascularization of the posterior coronary arteries had the following distribution: diagonalis artery (31 patients, 10%); marginal branches of the circumflex artery (243 patients, 78.7%); posterior ventricular artery (4 patients, 1.3%); and posterior descending artery (31 patients, 10%). RESULTS: Procedure-related complications without death occurred in 7 cases, giving a morbidity of 2.5%. There were 11 deaths in the early postoperative period (mortality of 3.9%). CONCLUSION: Similarly to the anterior coronary arteries, the posterior coronary arteries may benefit from myocardial revascularization without CPB.

Distribution of cardiac geometric patterns on echocardiography in essential hypertension. Impact of two criteria of stratification

PURPOSE: To evaluate 2 left ventricular mass index (LVMI) normality criteria for the prevalence of left ventricular geometric patterns in a hypertensive population ( HT ) . METHODS: 544 essential hypertensive patients, were evaluated by echocardiography, and different left ventricular hypertrophy criteria were applied: 1 - classic : men - 134 g/m² and women - 110 g/m² ; 2- obtained from the 95th percentil of LVMI from a normotensive population (NT). RESULTS: The prevalence of 4 left ventricular geometric patterns, respectively for criteria 1 and 2, were: normal geometry - 47.7% and 39.3%; concentric remodelying - 25.4% and 14.3%; concentric hypertrophy - 18.4% and 27.7% and excentric hypertrophy - 8.8% and 16.7%, which confered abnormal geometry to 52.6% and 60.7% of hypertensive. The comparative analysis between NT and normal geometry hypertensive group according to criteria 1, detected significative stuctural differences,"( *p < 0.05):LVMI- 78.4 ± 1.50 vs 85.9 ±0.95 g/m² *; posterior wall thickness -8.5 ± 0.1 vs 8.9 ± 0.05 mm*; left atrium - 33.3 ± 0.41 vs 34.7 ± 0.30 mm *. With criteria 2, significative structural differences between the 2 groups were not observed. CONCLUSION: The use of a reference population based criteria, increased the abnormal left ventricular geometry prevalence in hypertensive patients and seemed more appropriate for left ventricular hypertrophy detection and risk stratification.

Pharmacokinetics and posterior segment biodistribution of ESBA105, an anti-TNF-alpha single-chain antibody, upon topical administration to the rabbit eye.

PURPOSE: The purpose of this study was to characterize local distribution and systemic absorption of the tumor necrosis factor (TNF)-alpha inhibitory single-chain antibody fragment (scFv) ESBA105 following topical administration to the eye in vivo. METHODS: Rabbits received ESBA105 as topical eye drops in two dosing regimens. First, pharmacokinetics after the topical route of administration was compared to the intravenous (i.v.) route by means of applying the identical cumulative daily dose of ESBA105. In a second study rabbits received five eye drops daily for six consecutive days in a lower frequency topical dosing regimen. Kinetics and biodistribution of ESBA105 in ocular tissues and fluids as well as in sera were determined in all animals. RESULTS: After topical administration to the eye, ESBA105 quickly reaches therapeutic concentrations in all ocular compartments. Systemic exposure after topical administration is 25,000-fold lower than exposure after i.v. injection of the identical cumulative daily dose. ESBA105 levels in vitreous humor and neuroretina are significantly higher on topical administration than after i.v. injection. Absolute and relative intraocular biodistribution of ESBA105 is different with topical and systemic delivery routes. Compared to its terminal half-life in circulation (7 hours), the vitreal half-life of ESBA105 is significantly enhanced (16-24 hours). CONCLUSIONS: On topical administration, ESBA105 is efficiently absorbed and distributed to all compartments of the eye, whereby systemic drug exposure is very low. Based on its unique intraocular biodistribution and pharmacokinetics and the absolute intraocular levels reached, topical ESBA105 appears highly attractive for treatment of various ophthalmological disorders.