Showing papers in "Statistics and Computing in 2005"
TL;DR: This article describes how the series expansions can be summed in an numerically efficient fashion and demonstrates the usefulness of the approach, but full machine accuracy is shown not to be obtainable using the series expansion method for all parameter values.
Abstract: Exponential dispersion models, which are linear exponential families with a dispersion parameter, are the prototype response distributions for generalized linear models. The Tweedie family comprises those exponential dispersion models with power mean-variance relationships. The normal, Poisson, gamma and inverse Gaussian distributions belong to theTweedie family. Apart from these special cases, Tweedie distributions do not have density functions which can be written in closed form. Instead, the densities can be represented as infinite summations derived from series expansions. This article describes how the series expansions can be summed in an numerically efficient fashion. The usefulness of the approach is demonstrated, but full machine accuracy is shown not to be obtainable using the series expansion method for all parameter values. Derivatives of the density with respect to the dispersion parameter are also derived to facilitate maximum likelihood estimation. The methods are demonstrated on two data examples and compared with with Box-Cox transformations and extended quasi-likelihoood.
280 citations
TL;DR: In order to enlarge the applicability of the model, inference for a multivariate Poisson model with larger structure is proposed, i.e. different covariance for each pair of variables, and extension to models with complete structure with many multi-way covariance terms is discussed.
Abstract: In recent years the applications of multivariate Poisson models have increased, mainly because of the gradual increase in computer performance. The multivariate Poisson model used in practice is based on a common covariance term for all the pairs of variables. This is rather restrictive and does not allow for modelling the covariance structure of the data in a flexible way. In this paper we propose inference for a multivariate Poisson model with larger structure, i.e. different covariance for each pair of variables. Maximum likelihood estimation, as well as Bayesian estimation methods are proposed. Both are based on a data augmentation scheme that reflects the multivariate reduction derivation of the joint probability function. In order to enlarge the applicability of the model we allow for covariates in the specification of both the mean and the covariance parameters. Extension to models with complete structure with many multi-way covariance terms is discussed. The method is demonstrated by analyzing a real life data set.
143 citations
TL;DR: A Gaussian process mixture model for regression is proposed for dealing with the systematic heterogeneity among the different replications, and a hybrid Markov chain Monte Carlo (MCMC) algorithm is used for its implementation.
Abstract: As a result of their good performance in practice and their desirable analytical properties, Gaussian process regression models are becoming increasingly of interest in statistics, engineering and other fields. However, two major problems arise when the model is applied to a large data-set with repeated measurements. One stems from the systematic heterogeneity among the different replications, and the other is the requirement to invert a covariance matrix which is involved in the implementation of the model. The dimension of this matrix equals the sample size of the training data-set. In this paper, a Gaussian process mixture model for regression is proposed for dealing with the above two problems, and a hybrid Markov chain Monte Carlo (MCMC) algorithm is used for its implementation. Application to a real data-set is reported.
118 citations
TL;DR: Local principal curves are introduced, which are based on the localization of principal component analysis, and the proposed algorithm is able to identify closed curves as well as multiple curves which may or may not be connected.
Abstract: Principal components are a well established tool in dimension reduction. The extension to principal curves allows for general smooth curves which pass through the middle of a multidimensional data cloud. In this paper local principal curves are introduced, which are based on the localization of principal component analysis. The proposed algorithm is able to identify closed curves as well as multiple curves which may or may not be connected. For the evaluation of the performance of principal curves as tool for data reduction a measure of coverage is suggested. By use of simulated and real data sets the approach is compared to various alternative concepts of principal curves.
95 citations
TL;DR: This article applies Bayesian neural networks to time series analysis, and proposes a Monte Carlo algorithm for BNN training, and goes a step further in BNN model selection by putting a prior on network connections instead of hidden units as done by other authors.
Abstract: In this article, we apply Bayesian neural networks (BNNs) to time series analysis, and propose a Monte Carlo algorithm for BNN training In addition, we go a step further in BNN model selection by putting a prior on network connections instead of hidden units as done by other authors This allows us to treat the selection of hidden units and the selection of input variables uniformly The BNN model is compared to a number of competitors, such as the Box-Jenkins model, bilinear model, threshold autoregressive model, and traditional neural network model, on a number of popular and challenging data sets Numerical results show that the BNN model has achieved a consistent improvement over the competitors in forecasting future values Insights on how to improve the generalization ability of BNNs are revealed in many respects of our implementation, such as the selection of input variables, the specification of prior distributions, and the treatment of outliers
94 citations
TL;DR: It is shown through simulation that no single model selection criterion exhibits a uniformly superior performance over a wide range of scenarios, so a two-stage approach for model selection is proposed and shown to perform satisfactorily.
Abstract: The performance of computationally inexpensive model selection criteria in the context of tree-structured subgroup analysis is investigated. It is shown through simulation that no single model selection criterion exhibits a uniformly superior performance over a wide range of scenarios. Therefore, a two-stage approach for model selection is proposed and shown to perform satisfactorily. Applied example of subgroup analysis is presented. Problems associated with tree-structured subgroup analysis are discussed and practical solutions are suggested.
72 citations
TL;DR: Non-centered and partially non-centered MCMC algorithms for stochastic epidemic models are introduced and are shown to out perform the existing centered algorithms.
Abstract: In this paper, we introduce non-centered and partially non-centered MCMC algorithms for stochastic epidemic models. Centered algorithms previously considered in the literature perform adequately well for small data sets. However, due to the high dependence inherent in the models between the missing data and the parameters, the performance of the centered algorithms gets appreciably worse when larger data sets are considered. Therefore non-centered and partially non-centered algorithms are introduced and are shown to out perform the existing centered algorithms.
71 citations
TL;DR: Several new sequential Monte Carlo algorithms for online estimation (filtering) of nonlinear dynamic systems and are efficient because they tend to utilize both the information in the state process and the observations and are easy to sample from.
Abstract: In this paper we present several new sequential Monte Carlo (SMC) algorithms for online estimation (filtering) of nonlinear dynamic systems. SMC has been shown to be a powerful tool for dealing with complex dynamic systems. It sequentially generates Monte Carlo samples from a proposal distribution, adjusted by a set of importance weight with respect to a target distribution, to facilitate statistical inferences on the characteristic (state) of the system. The key to a successful implementation of SMC in complex problems is the design of an efficient proposal distribution from which the Monte Carlo samples are generated. We propose several such proposal distributions that are efficient yet easy to generate samples from. They are efficient because they tend to utilize both the information in the state process and the observations. They are all Gaussian distributions hence are easy to sample from. The central ideas of the conventional nonlinear filters, such as extended Kalman filter, unscented Kalman filter and the Gaussian quadrature filter, are used to construct these proposal distributions. The effectiveness of the proposed algorithms are demonstrated through two applications--real time target tracking and the multiuser parameter tracking in CDMA communication systems.
61 citations
TL;DR: This work exploits the collocation method and the product Nyström method and sees that collocation leads to higher accuracy than currently established methods.
Abstract: Originally, the exponentially weighted moving average (EWMA) control chart was developed for detecting changes in the process mean. The average run length (ARL) became the most popular performance measure for schemes with this objective. When monitoring the mean of independent and normally distributed observations the ARL can be determined with high precision. Nowadays, EWMA control charts are also used for monitoring the variance. Charts based on the sample variance S2 are an appropriate choice. The usage of ARL evaluation techniques known from mean monitoring charts, however, is difficult. The most accurate method--solving a Fredholm integral equation with the Nystrom method--fails due to an improper kernel in the case of chi-squared distributions. Here, we exploit the collocation method and the product Nystrom method. These methods are compared to Markov chain based approaches. We see that collocation leads to higher accuracy than currently established methods.
57 citations
TL;DR: An auxiliary variable method based on a slice sampler is shown to provide an attractive simulation-based model fitting strategy for fitting Bayesian models under proper priors.
Abstract: An auxiliary variable method based on a slice sampler is shown to provide an attractive simulation-based model fitting strategy for fitting Bayesian models under proper priors. Though broadly applicable, we illustrate in the context of fitting spatial models for geo-referenced or point source data. Spatial modeling within a Bayesian framework offers inferential advantages and the slice sampler provides an algorithm which is essentially "off the shelf". Further potential advantages over importance sampling approaches and Metropolis approaches are noted and illustrative examples are supplied.
52 citations
TL;DR: Optimal bandwidth choice for matching estimators and their finite sample properties are examined and an approximation to their MSE is derived, as a basis for a plug-in bandwidth selector.
Abstract: Optimal bandwidth choice for matching estimators and their finite sample properties are examined. An approximation to their MSE is derived, as a basis for a plug-in bandwidth selector. In small samples, this approximation is not very accurate, though. Alternatively, conventional cross-validation bandwidth selection is considered and performs rather well in simulation studies: Compared to standard pair-matching, kernel and ridge matching achieve reductions in MSE of about 25 to 40%. Local linear matching and weighting perform poorly. Furthermore, the scope for developing better bandwidth selectors seems to be limited for ridge matching, but non-negligible for kernel and local linear matching.
TL;DR: Connections between the frequentist P-value and the posterior distribution of the likelihood ratio are used to interpret and calibrate P-values in a Bayesian context, and examples are given to show the use of simple posterior simulation methods to provide Bayesian tests of common hypotheses.
Abstract: This paper gives an exposition of the use of the posterior likelihood ratio for testing point null hypotheses in a fully Bayesian framework. Connections between the frequentist P-value and the posterior distribution of the likelihood ratio are used to interpret and calibrate P-values in a Bayesian context, and examples are given to show the use of simple posterior simulation methods to provide Bayesian tests of common hypotheses.
TL;DR: It is shown that boosting kernel classifiers reduces the bias whilst only slightly increasing the variance, with an overall reduction in error, which is closely linked to a previously proposed method of bias reduction in kernel density estimation.
Abstract: Kernel density estimation is a commonly used approach to classification. However, most of the theoretical results for kernel methods apply to estimation per se and not necessarily to classification. In this paper we show that when estimating the difference between two densities, the optimal smoothing parameters are increasing functions of the sample size of the complementary group, and we provide a small simluation study which examines the relative performance of kernel density methods when the final goal is classification.
A relative newcomer to the classification portfolio is "boosting", and this paper proposes an algorithm for boosting kernel density classifiers. We note that boosting is closely linked to a previously proposed method of bias reduction in kernel density estimation and indicate how it will enjoy similar properties for classification. We show that boosting kernel classifiers reduces the bias whilst only slightly increasing the variance, with an overall reduction in error. Numerical examples and simulations are used to illustrate the findings, and we also suggest further areas of research.
TL;DR: This article proposes a fast recursive construction of the inner product matrix of discrete a.c. wavelets which is required by the statistical analysis and describes an efficient construction in the (separable) two-dimensional case.
Abstract: Discrete autocorrelation (a.c.) wavelets have recently been applied in the statistical analysis of locally stationary time series for local spectral modelling and estimation. This article proposes a fast recursive construction of the inner product matrix of discrete a.c. wavelets which is required by the statistical analysis. The recursion connects neighbouring elements on diagonals of the inner product matrix using a two-scale property of the a.c. wavelets. The recursive method is an ?(log (N)3) operation which compares favourably with the ?(N log N) operations required by the brute force approach. We conclude by describing an efficient construction of the inner product matrix in the (separable) two-dimensional case.
TL;DR: A local scoring algorithm (with backfitting) based on local linear kernel smoothers was used to estimate a generalized additive model with second-order interaction terms and a bootstrap procedure is provided for estimating the distribution of the test statistics.
Abstract: In this paper we considered a generalized additive model with second-order interaction terms. A local scoring algorithm (with backfitting) based on local linear kernel smoothers was used to estimate the model. Our main aim was to obtain procedures for testing second-order interaction terms. Backfitting theory is difficult in this context, and a bootstrap procedure is therefore provided for estimating the distribution of the test statistics. Given the high computational cost involved, binning techniques were used to speed up the computation in the estimation and testing process. A simulation study was carried out in order to assess the validity of the bootstrap-based tests. Lastly, our method was applied to real data drawn from an SO2 binary time series.
TL;DR: The Unobserved ARCH model is a good description of the phenomenon of changing volatility that is commonly appeared in the financial time series and some suitable non-linear transformations of the parameter space are adopted such that the resulting MCMC algorithm is based only on Gibbs sampling steps.
Abstract: The Unobserved ARCH model is a good description of the phenomenon of changing volatility that is commonly appeared in the financial time series We study this model adopting Bayesian inference via Markov Chain Monte Carlo (MCMC) In order to provide an easy to implement MCMC algorithm we adopt some suitable non-linear transformations of the parameter space such that the resulting MCMC algorithm is based only on Gibbs sampling steps We illustrate our methodology with data from real world The Unobserved ARCH is shown to be a good description of the exchange rate movements Numerical comparisons between competing MCMC algorithms are also presented
TL;DR: The approach is based on directly calculating the posterior distribution using a set of recursions which are similar to those of the Forward-Backward algorithm, which is more practicable than existing perfect simulation methods for mixtures.
Abstract: We demonstrate how to perform direct simulation for discrete mixture models. The approach is based on directly calculating the posterior distribution using a set of recursions which are similar to those of the Forward-Backward algorithm. Our approach is more practicable than existing perfect simulation methods for mixtures. For example, we analyse 1096 observations from a 2 component Poisson mixture, and 240 observations under a 3 component Poisson mixture (with unknown mixture proportions and Poisson means in each case). Simulating samples of 10,000 perfect realisations took about 17 minutes and an hour respectively on a 900 MHz ultraSPARC computer. Our method can also be used to perform perfect simulation from Markov-dependent mixture models. A byproduct of our approach is that the evidence of our assumed models can be calculated, which enables different models to be compared.
TL;DR: The choice of the GEP density as an importance function allows us to obtain reliable and effective results when p-credences of the prior and the likelihood are defined, even if there are conflicting sources of information.
Abstract: In this paper, the generalized exponential power (GEP) density is proposed as an importance function in Monte Carlo simulations in the context of estimation of posterior moments of a location parameter. This density is divided in five classes according to its tail behaviour which may be exponential, polynomial or logarithmic. The notion of p-credence is also defined to characterize and to order the tails of a large class of symmetric densities by comparing their tails to those of the GEP density.
The choice of the GEP density as an importance function allows us to obtain reliable and effective results when p-credences of the prior and the likelihood are defined, even if there are conflicting sources of information. Characterization of the posterior tails using p-credence can be done. Hence, it is possible to choose parameters of the GEP density in order to have an importance function with slightly heavier tails than the posterior. Simulation of observations from the GEP density is also addressed.
TL;DR: This paper suggests a similar idea in which the Metropolis-Hastings proposals of Denison, Mallick and Smith (1998a) are altered to allow dependence on the current model, which allows more rapid identification and exploration of important interactions, especially in problems with very large numbers of predictor variables and many useless predictors.
Abstract: Multivariate adaptive regression spline fitting or MARS (Friedman 1991) provides a useful methodology for flexible adaptive regression with many predictors. The MARS methodology produces an estimate of the mean response that is a linear combination of adaptively chosen basis functions. Recently, a Bayesian version of MARS has been proposed (Denison, Mallick and Smith 1998a, Holmes and Denison, 2002) combining the MARS methodology with the benefits of Bayesian methods for accounting for model uncertainty to achieve improvements in predictive performance. In implementation of the Bayesian MARS approach, Markov chain Monte Carlo methods are used for computations, in which at each iteration of the algorithm it is proposed to change the current model by either (a) Adding a basis function (birth step) (b) Deleting a basis function (death step) or (c) Altering an existing basis function (change step). In the algorithm of Denison, Mallick and Smith (1998a), when a birth step is proposed, the type of basis function is determined by simulation from the prior. This works well in problems with a small number of predictors, is simple to program, and leads to a simple form for Metropolis-Hastings acceptance probabilities. However, in problems with very large numbers of predictors where many of the predictors are useless it may be difficult to find interesting interactions with such an approach. In the original MARS algorithm of Friedman (1991) a heuristic is used of building up higher order interactions from lower order ones, which greatly reduces the complexity of the search for good basis functions to add to the model. While we do not exactly follow the intuition of the original MARS algorithm in this paper, we nevertheless suggest a similar idea in which the Metropolis-Hastings proposals of Denison, Mallick and Smith (1998a) are altered to allow dependence on the current model. Our modification allows more rapid identification and exploration of important interactions, especially in problems with very large numbers of predictor variables and many useless predictors. Performance of the algorithms is compared in simulation studies.
TL;DR: This paper presents simple yet extremely accurate saddlepoint approximations to power functions associated with the following classical test statistics: the likelihood ratio statistic for testing the general linear hypothesis in MANOVA; the likelihood ratios for testing block independence; and Bartlett's modified likelihood ratio statistics for testing equality of covariance matrices.
Abstract: We consider the calculation of power functions in classical multivariate analysis. In this context, power can be expressed in terms of tail probabilities of certain noncentral distributions. The necessary noncentral distribution theory was developed between the 1940s and 1970s by a number of authors. However, tractable methods for calculating the relevant probabilities have been lacking. In this paper we present simple yet extremely accurate saddlepoint approximations to power functions associated with the following classical test statistics: the likelihood ratio statistic for testing the general linear hypothesis in MANOVA; the likelihood ratio statistic for testing block independence; and Bartlett's modified likelihood ratio statistic for testing equality of covariance matrices.
TL;DR: The number of components k is obtained as the mode of the bootstrap distribution of k, which is presented using the Times newspaper data and investigated in a simulation study for mixtures of Poisson data.
Abstract: Finite mixture models arise in a natural way in that they are modeling unobserved population heterogeneity. It is assumed that the population consists of an unknown number k of subpopulations with parameters ?1, ..., ?k receiving weights p1, ..., pk. Because of the irregularity of the parameter space, the log-likelihood-ratio statistic (LRS) does not have a (?2) limit distribution and therefore it is difficult to use the LRS to test for the number of components. These problems are circumvented by using the nonparametric bootstrap such that the mixture algorithm is applied B times to bootstrap samples obtained from the original sample with replacement. The number of components k is obtained as the mode of the bootstrap distribution of k. This approach is presented using the Times newspaper data and investigated in a simulation study for mixtures of Poisson data.
TL;DR: It is shown that the P-value resulting from the hypothesis test, considered as a function of the null-hypothesized value of θ, has both “jump” and “drop” discontinuities.
Abstract: A new area of research interest is the computation of exact confidence limits or intervals for a scalar parameter of interest ? from discrete data by inverting a hypothesis test based on a studentized test statistic. See, for example, Chan and Zhang (1999), Agresti and Min (2001) and Agresti (2003) who deal with ? a difference of binomial probabilities and Agresti and Min (2002) who deal with ? an odds ratio. However, neither (1) a detailed analysis of the computational issues involved nor (2) a reliable method of computation that deals effectively with these issues is currently available. In this paper we solve these two problems for a very broad class of discrete data models. We suppose that the distribution of the data is determined by (?,?) where ? is a nuisance parameter vector. We also consider six different studentized test statistics. Our contributions to (1) are as follows. We show that the P-value resulting from the hypothesis test, considered as a function of the null-hypothesized value of ?, has both "jump" and "drop" discontinuities. Numerical examples are used to demonstrate that these discontinuities lead to the failure of simple-minded approaches to the computation of the confidence limit or interval. We also provide a new method for efficiently computing the set of all possible locations of these discontinuities. Our contribution to (2) is to provide a new and reliable method of computing the confidence limit or interval, based on the knowledge of this set.
TL;DR: It is shown that this posterior predictive distribution formula derived in Sweeting And Kharroubi (2003) provides a stable importance function for use within poor man’s data augmentation schemes and that it can also be used as a proposal distribution within a Metropolis-Hastings algorithm for models that are not analytically tractable.
Abstract: We consider exact and approximate Bayesian computation in the presence of latent variables or missing data. Specifically we explore the application of a posterior predictive distribution formula derived in Sweeting And Kharroubi (2003), which is a particular form of Laplace approximation, both as an importance function and a proposal distribution. We show that this formula provides a stable importance function for use within poor man's data augmentation schemes and that it can also be used as a proposal distribution within a Metropolis-Hastings algorithm for models that are not analytically tractable. We illustrate both uses in the case of a censored regression model and a normal hierarchical model, with both normal and Student t distributed random effects. Although the predictive distribution formula is motivated by regular asymptotic theory, it is not necessary that the likelihood has a closed form or that it possesses a local maximum.
TL;DR: The task of determining expected values of sample moments, where the sample members have been selected based on noisy information, is considered and it is shown experimentally that including skewness and kurtosis in the calculations can yield greatly improved results for other distributions.
Abstract: In this paper, the task of determining expected values of sample moments, where the sample members have been selected based on noisy information, is considered. This task is a recurring problem in the theory of evolution strategies. Exact expressions for expected values of sums of products of concomitants of selected order statistics are derived. Then, using Edgeworth and Cornish-Fisher approximations, explicit results that depend on coefficients that can be determined numerically are obtained. While the results are exact only for normal populations, it is shown experimentally that including skewness and kurtosis in the calculations can yield greatly improved results for other distributions.
TL;DR: An important merit of the proposed Yao-Tong approach is that it is conceptually simple and can be readily applied to parametrically nonlinear conditional quantile estimation.
Abstract: In this paper, nonparametric estimation of conditional quantiles of a nonlinear time series model is formulated as a nonsmooth optimization problem involving an asymmetric loss function. This asymmetric loss function is nonsmooth and is of the same structure as the so-called `lopsided' absolute value function. Using an effective smoothing approximation method introduced for this lopsided absolute value function, we obtain a sequence of approximate smooth optimization problems. Some important convergence properties of the approximation are established. Each of these smooth approximate optimization problems is solved by an optimization algorithm based on a sequential quadratic programming approximation with active set strategy. Within the framework of locally linear conditional quantiles, the proposed approach is compared with three other approaches, namely, an approach proposed by Yao and Tong (1996), the Iteratively Reweighted Least Squares method and the Interior-Point method, through some empirical numerical studies using simulated data and the classic lynx pelt series. In particular, the empirical performance of the proposed approach is almost identical with that of the Interior-Point method, both methods being slightly better than the Iteratively Reweighted Least Squares method. The Yao-Tong approach is comparable with the other methods in the ideal cases for the Yao-Tong method, but otherwise it is outperformed by other approaches. An important merit of the proposed approach is that it is conceptually simple and can be readily applied to parametrically nonlinear conditional quantile estimation.
TL;DR: Experimental studies demonstrate that this Bayesian source separation algorithm is appropriate for systematic spatial pattern analysis by modeling arbitrary sources and identify their effects on high dimensional measurement data.
Abstract: A Bayesian blind source separation (BSS) algorithm is proposed in this paper to recover independent sources from observed multivariate spatial patterns. As a widely used mechanism, Gaussian mixture model is adopted to represent the sources for statistical description and machine learning. In the context of linear latent variable BSS model, some conjugate priors are incorporated into the hyperparameters estimation of mixing matrix. The proposed algorithm then approximates the full posteriors over model structure and source parameters in an analytical manner based on variational Bayesian treatment. Experimental studies demonstrate that this Bayesian source separation algorithm is appropriate for systematic spatial pattern analysis by modeling arbitrary sources and identify their effects on high dimensional measurement data. The identified patterns will serve as diagnosis aids for gaining insight into the nature of physical process for the potential use of statistical quality control.
TL;DR: This paper first develops relevant dynamical and communication theory in the bivariate map context, and then presents a better way of improving synchronization by distribution transformation, which involves transforming the transmission sequence using knowledge of the invariant distribution of the spreading sequence, and before noise corrupts the signal in the transmission channel.
Abstract: Research in electronic communications has developed chaos-based modelling to enable messages to be carried by chaotic broad-band spreading sequences. When such systems are used it is necessary to simultaneously know the spreading sequence at both the transmitting and receiving stations. This is possible using the idea of synchronization with bivariate maps, providing there is no noise present in the system. When noise is present in the transmission channel, recovery of the spreading sequence may be degraded or impossible. Once noise is added to the spreading sequence, the result may no longer lie within the boundary of the chaotic map. A usual and obvious method of dealing with this problem is to cap iterations lying outside the bounds at their extremes, but the procedure amplifies loss of synchronization. With a minimum of technical details and a computational focus, this paper first develops relevant dynamical and communication theory in the bivariate map context, and then presents a better way of improving synchronization by distribution transformation. The transmission sequence is transformed, using knowledge of the invariant distribution of the spreading sequence, and before noise corrupts the signal in the transmission channel. An `inverse' transformation can then be applied at the receiver station so that the noise has a reduced impact on the recovery of the spreading sequence and hence its synchronization. Statistical simulations illustrating the effectiveness of the approach are presented.