Likelihood Surface Research Articles

Neural likelihood surfaces for spatial processes with computationally intensive or intractable likelihoods

In spatial statistics, fast and accurate parameter estimation, coupled with a reliable means of uncertainty quantification, can be challenging when fitting a spatial process to real-world data because the likelihood function might be slow to evaluate or wholly intractable. In this work, we propose using convolutional neural networks to learn the likelihood function of a spatial process. Through a specifically designed classification task, our neural network implicitly learns the likelihood function, even in situations where the exact likelihood is not explicitly available. Once trained on the classification task, our neural network is calibrated using Platt scaling which improves the accuracy of the neural likelihood surfaces. To demonstrate our approach, we compare neural likelihood surfaces and the resulting maximum likelihood estimates and approximate confidence regions with the equivalent for exact or approximate likelihood for two different spatial processes—a Gaussian process and a Brown–Resnick process which have computationally intensive and intractable likelihoods, respectively. We conclude that our method provides fast and accurate parameter estimation with a reliable method of uncertainty quantification in situations where standard methods are either undesirably slow or inaccurate. The method is applicable to any spatial process on a grid from which fast simulations are available.

Estimating Galaxy Parameters with Self-organizing Maps and the Effect of Missing Data

The current and upcoming large data volume galaxy surveys require the use of machine-learning techniques to maximize their scientific return. This study explores the use of Self-Organizing Maps (SOMs) to estimate galaxy parameters with a focus on handling cases of missing data and providing realistic probability distribution functions for the parameters. We train an SOM with a simulated mass-limited lightcone assuming a ugrizY JHK s +IRAC data set, mimicking the Hyper Suprime-Cam Deep joint data set. For parameter estimation, we derive SOM likelihood surfaces considering photometric errors to derive total (statistical and systematic) uncertainties. We explore the effects of missing data, including which bands are particularly critical to the accuracy of the derived parameters. We demonstrate that the parameter recovery is significantly better when the missing bands are “filled in” rather than if they are completely omitted. We propose a practical method for such recovery of missing data.

Understanding the impact of numerical solvers on inference for differential equation models.

Most ordinary differential equation (ODE) models used to describe biological or physical systems must be solved approximately using numerical methods. Perniciously, even those solvers that seem sufficiently accurate for the forward problem, i.e. for obtaining an accurate simulation, might not be sufficiently accurate for the inverse problem, i.e. for inferring the model parameters from data. We show that for both fixed step and adaptive step ODE solvers, solving the forward problem with insufficient accuracy can distort likelihood surfaces, which might become jagged, causing inference algorithms to get stuck in local 'phantom' optima. We demonstrate that biases in inference arising from numerical approximation of ODEs are potentially most severe in systems involving low noise and rapid nonlinear dynamics. We reanalyse an ODE changepoint model previously fit to the COVID-19 outbreak in Germany and show the effect of the step size on simulation and inference results. We then fit a more complicated rainfallrun-off model to hydrological data and illustrate the importance of tuning solver tolerances to avoid distorted likelihood surfaces. Our results indicate that, when performing inference for ODE model parameters, adaptive step size solver tolerances must be set cautiously and likelihood surfaces should be inspected for characteristic signs of numerical issues.

Dentist: Quantifying uncertainty by sampling points around maximum likelihood estimates

Abstract It is standard statistical practice to provide measures of uncertainty around parameter estimates. Unfortunately, this very basic and necessary enterprise is often absent in macroevolutionary studies using maximum likelihood estimates (MLEs). dentist is an R package that allows an approximation of confidence intervals (CI) around parameter estimates without an analytic solution to likelihood equations. This package works by ‘denting’ the likelihood surface by sampling points a specified distance around the MLE following what is essentially a Metropolis‐Hastings walk. We describe the importance of estimating uncertainty around parameter estimates, as well as demonstrate the ability of dentist to accurately approximate CI. We introduce several plotting tools to visualize the results of a dentist analysis. dentist is freely available from https://github.com/bomeara/dentist, written in the R language, and can be used for any given likelihood function.

Adaptive RAxML-NG: Accelerating Phylogenetic Inference under Maximum Likelihood using Dataset Difficulty.

Phylogenetic inferences under the maximum likelihood criterion deploy heuristic tree search strategies to explore the vast search space. Depending on the input dataset, searches from different starting trees might all converge to a single tree topology. Often, though, distinct searches infer multiple topologies with large log-likelihood score differences or yield topologically highly distinct, yet almost equally likely, trees. Recently, Haag et al. introduced an approach to quantify, and implemented machine learning methods to predict, the dataset difficulty with respect to phylogenetic inference. Easy multiple sequence alignments (MSAs) exhibit a single likelihood peak on their likelihood surface, associated with a single tree topology to which most, if not all, independent searches rapidly converge. As difficulty increases, multiple locally optimal likelihood peaks emerge, yet from highly distinct topologies. To make use of this information, we introduce and implement an adaptive tree search heuristic in RAxML-NG, which modifies the thoroughness of the tree search strategy as a function of the predicted difficulty. Our adaptive strategy is based upon three observations. First, on easy datasets, searches converge rapidly and can hence be terminated at an earlier stage. Second, overanalyzing difficult datasets is hopeless, and thus it suffices to quickly infer only one of the numerous almost equally likely topologies to reduce overall execution time. Third, more extensive searches are justified and required on datasets with intermediate difficulty. While the likelihood surface exhibits multiple locally optimal peaks in this case, a small proportion of them is significantly better. Our experimental results for the adaptive heuristic on 9,515 empirical and 5,000 simulated datasets with varying difficulty exhibit substantial speedups, especially on easy and difficult datasets (53% of total MSAs), where we observe average speedups of more than 10×. Further, approximately 94% of the inferred trees using the adaptive strategy are statistically indistinguishable from the trees inferred under the standard strategy (RAxML-NG).

On near-redundancy and identifiability of parametric hazard regression models under censoring.

We study parametric inference on a rich class of hazard regression models in the presence of right-censoring. Previous literature has reported some inferential challenges, such as multimodal or flat likelihood surfaces, in this class of models for some particular data sets. We formalize the study of these inferential problems by linking them to the concepts of near-redundancy and practical nonidentifiability of parameters. We show that the maximum likelihood estimators of the parameters in this class of models are consistent and asymptotically normal. Thus, the inferential problems in this class of models are related to the finite-sample scenario, where it is difficult to distinguish between the fitted model and a nested nonidentifiable (i.e., parameter-redundant) model. We propose a method for detecting near-redundancy, based on distances between probability distributions. We also employ methods used in other areas for detecting practical nonidentifiability and near-redundancy, including the inspection of the profile likelihood function and the Hessian method. For cases where inferential problems are detected, we discuss alternatives such as using model selection tools to identify simpler models that do not exhibit these inferential problems, increasing the sample size, or extending the follow-up time. We illustrate the performance of the proposed methods through a simulation study. Our simulation study reveals a link between the presence of near-redundancy and practical nonidentifiability. Two illustrative applications using real data, with and without inferential problems, arepresented.

Halo occupation distribution of Emission Line Galaxies: fitting method with Gaussian processes

The halo occupation distribution (HOD) framework is an empirical method to describe the connection between dark matter halos and galaxies, which is constrained by small scale clustering data. Efficient fitting procedures are required to scan the HOD parameter space. This paper describes such a method based on Gaussian Processes to iteratively build a surrogate model of the posterior of the likelihood surface from a reasonable amount of likelihood computations, typically two orders of magnitude less than standard Monte Carlo Markov chain algorithms. Errors in the likelihood computation due to stochastic HOD modelling are also accounted for in the method we propose. We report results of reproducibility, accuracy and stability tests of the method derived from simulation, taking as a test case star-forming emission line galaxies, which constitute the main tracer of the Dark Energy Spectroscopic Instrument and have so far a poorly constrained galaxy-halo connection from observational data.

On the parameter estimation of Box-Cox transformation cure model.

We propose an improved estimation method for the Box-Cox transformation (BCT) cure rate model parameters. Specifically, we propose a generic maximum likelihood estimation algorithm through a non-linear conjugate gradient (NCG) method with an efficient line search technique. We then apply the proposed NCG algorithm to BCT cure model. Through a detailed simulation study, we compare the model fitting results of the NCG algorithm with those obtained by the existing expectation maximization (EM) algorithm. First, we show that our proposed NCG algorithm allows simultaneous maximization of all model parameters unlike the EM algorithm when the likelihood surface is flat with respect to the BCT index parameter. Then, we show that the NCG algorithm results in smaller bias and noticeably smaller root mean square error of the estimates of the model parameters that are associated with the cure rate. This results in more accurate and precise inference on the cure rate. In addition, we show that when the sample size is large the NCG algorithm, which only needs the computation of the gradient and not the Hessian, takes less CPU time to produce the estimates. These advantages of the NCG algorithm allows us to conclude that the NCG method should be the preferred estimation method over the already existing EM algorithm in the context of BCT cure model. Finally, we apply the NCG algorithm to analyze a well-known melanoma data and show that it results in a better fit when compared to the EM algorithm.

Combined full shape analysis of BOSS galaxies and eBOSS quasars using an iterative emulator

ABSTRACT Standard full-shape clustering analyses in Fourier space rely on a fixed power spectrum template, defined at the fiducial cosmology used to convert redshifts into distances, and compress the cosmological information into the Alcock–Paczynski parameters and the linear growth rate of structure. In this paper, we propose an analysis method that operates directly in the cosmology parameter space and varies the power spectrum template accordingly at each tested point. Predictions for the power spectrum multipoles from the TNS model are computed at different cosmologies in the framework of $\Lambda \rm {CDM}$. Applied to the final eBOSS QSO and LRG samples together with the low-z DR12 BOSS galaxy sample, our analysis results in a set of constraints on the cosmological parameters Ωcdm, H0, σ8, Ωb, and ns. To reduce the number of computed models, we construct an iterative process to sample the likelihood surface, where each iteration consists of a Gaussian process regression. This method is validated with mocks from N-body simulations. From the combined analysis of the (e)BOSS data, we obtain the following constraints: σ8 = 0.877 ± 0.049 and $\Omega _{\rm m}=0.304^{+0.016}_{-0.010}$ without any external prior. The eBOSS quasar sample alone shows a 3.1σ discrepancy compared to the Planck prediction.

Density reconstruction in convergent high-energy-density systems using x-ray radiography and Bayesian inference

X-ray radiography is a technique frequently used to diagnose convergent high-energy-density (HED) systems, such as inertial confinement fusion implosion, and provides unique information that is not available through self-emission measurements. We investigate the scope and limits of that information using a radiography simulation combined with the Bayesian inference workflow. The accuracy of density reconstruction from simulated radiographs of spherical implosions driven with 27 kJ laser energy is assessed, including the increase or decrease in accuracy due to the addition of Lagrangian marker layers, Poisson noise, and improved prior information. This work is the first to present the full uncertainty distributions inferred from radiography analysis in HED systems and demonstrates the importance of constructing the full posterior probability density, as opposed to a point estimate, due to the modal structure of the likelihood surface introduced by typical experimental noise sources. This general methodology can be used both for robust analysis of radiographic data and for an improved design of radiography experiments by modeling the full experimental system.

Probabilistic seismic source location and magnitude via inverse analysis of paleoliquefaction evidence

In regions of infrequent moderate-to-large earthquakes, historic earthquake catalogs are often insufficient to provide inputs to seismic-hazard analyses (i.e. fault locations and magnitude–frequency relations) or to inform ground-motion predictions for certain seismic sources. In these regions, analysis of relic coseismic evidence, such as paleoliquefaction, is commonly used to infer information about the seismic hazard. However, while paleoliquefaction studies have been performed widely, all existing analysis techniques require a priori assumptions about the causative earthquake’s location (i.e. rupture magnitude and ground motions can otherwise not be estimated). This may lead to inaccurate assumptions in some settings, and by corollary, erroneous results. Accordingly, this article proposes an inversion framework to probabilistically constrain seismic-source parameters from paleoliquefaction. Analyzing evidence at regional scale leads to (a) a geospatial likelihood surface that constrains the rupture location and (b) a probability distribution of the rupture magnitude, wherein source-location uncertainty is explicitly considered. Simulated paleoliquefaction studies are performed on earthquakes with known parameters. These examples demonstrate the framework’s potential, even in cases of limited field evidence, as well as important caveats and lessons for forward use. The proposed framework has the potential to provide new insights in enigmatic seismic zones worldwide.

Adaptive Semiparametric Bayesian Differential Equations Via Sequential Monte Carlo

Nonlinear differential equations (DEs) are used in a wide range of scientific problems to model complex dynamic systems. The differential equations often contain unknown parameters that are of scientific interest, which have to be estimated from noisy measurements of the dynamic system. Generally, there is no closed-form solution for nonlinear DEs, and the likelihood surface for the parameter of interest is multi-modal and very sensitive to different parameter values. We propose a Bayesian framework for nonlinear DE systems. A flexible nonparametric function is used to represent the dynamic process such that expensive numerical solvers can be avoided. A sequential Monte Carlo algorithm in the annealing framework is proposed to conduct Bayesian inference for parameters in DEs. In our numerical experiments, we use examples of ordinary differential equations and delay differential equations to demonstrate the effectiveness of the proposed algorithm. We developed an R package that is available at https://github.com/shijiaw/smcDE. Supplementary files for this article are available online.

Estimation of Natural Selection and Allele Age from Time Series Allele Frequency Data Using a Novel Likelihood-Based Approach.

Temporally spaced genetic data allow for more accurate inference of population genetic parameters and hypothesis testing on the recent action of natural selection. In this work, we develop a novel likelihood-based method for jointly estimating selection coefficient and allele age from time series data of allele frequencies. Our approach is based on a hidden Markov model where the underlying process is a Wright-Fisher diffusion conditioned to survive until the time of the most recent sample. This formulation circumvents the assumption required in existing methods that the allele is created by mutation at a certain low frequency. We calculate the likelihood by numerically solving the resulting Kolmogorov backward equation backward in time while reweighting the solution with the emission probabilities of the observation at each sampling time point. This procedure reduces the two-dimensional numerical search for the maximum of the likelihood surface, for both the selection coefficient and the allele age, to a one-dimensional search over the selection coefficient only. We illustrate through extensive simulations that our method can produce accurate estimates of the selection coefficient and the allele age under both constant and nonconstant demographic histories. We apply our approach to reanalyze ancient DNA data associated with horse base coat colors. We find that ignoring demographic histories or grouping raw samples can significantly bias the inference results.

An Evolutionary Algorithm with Crossover and Mutation for Model-Based Clustering

An evolutionary algorithm (EA) is developed as an alternative to the EM algorithm for parameter estimation in model-based clustering. This EA facilitates a different search of the fitness landscape, i.e., the likelihood surface, utilizing both crossover and mutation. Furthermore, this EA represents an efficient approach to “hard” model-based clustering and so it can be viewed as a sort of generalization of the k-means algorithm, which is itself equivalent to a restricted Gaussian mixture model. The EA is illustrated on several datasets, and its performance is compared with that of other hard clustering approaches and model-based clustering via the EM algorithm.

Penalized variance components for association of multiple genes with traits.

Variance component models have gained popularity for genetic analyses, driven by their flexibility to simultaneously analyze multiple genetic variants in a gene by kernel statistics, and their ability to account for population stratification via genomic relationship matrices. For exploratory analyses with modest sample sizes and a potentially large number of variance components, it can be challenging to use standard maximum-likelihood or restricted maximum-likelihood methods to estimate variance components, because these iterative methods often fail to converge when likelihood surfaces are fairly flat, and standard-likelihood ratio statistical tests are not adequate. To overcome these limitations, we developed a penalized-likelihood model, whereby the penalty function follows the popular elastic-net approach, applying both L1 and L2 penalties to the variance components. By simulations, we demonstrate the potential gain in power by using both L1 and L2 penalties, and results from our simulations suggest that assigning 80% of the penalty parameter to the L1 penalty and 20% to the L2 penalty provides a reasonable balance between false-positive and false-negative results. Larger sample size improves the properties of our methods, at the cost of longer computation time. Application of our methods to a study of the influence of DNA methylation on levels of cortisol in reaction to stress testing shows how our method can be used to prioritize findings for further functional studies.

Saddle-Reset for Robust Parameter Estimation and Identifiability Analysis of Nonlinear Mixed Effects Models

Parameter estimation of a nonlinear model based on maximizing the likelihood using gradient-based numerical optimization methods can often fail due to premature termination of the optimization algorithm. One reason for such failure is that these numerical optimization methods cannot distinguish between the minimum, maximum, and a saddle point; hence, the parameters found by these optimization algorithms can possibly be in any of these three stationary points on the likelihood surface. We have found that for maximization of the likelihood for nonlinear mixed effects models used in pharmaceutical development, the optimization algorithm Broyden–Fletcher–Goldfarb–Shanno (BFGS) often terminates in saddle points, and we propose an algorithm, saddle-reset, to avoid the termination at saddle points, based on the second partial derivative test. In this algorithm, we use the approximated Hessian matrix at the point where BFGS terminates, perturb the point in the direction of the eigenvector associated with the lowest eigenvalue, and restart the BFGS algorithm. We have implemented this algorithm in industry standard software for nonlinear mixed effects modeling (NONMEM, version 7.4 and up) and showed that it can be used to avoid termination of parameter estimation at saddle points, as well as unveil practical parameter non-identifiability. We demonstrate this using four published pharmacometric models and two models specifically designed to be practically non-identifiable.

Epigenetic pacemaker: closed form algebraic solutions

BackgroundDNA methylation is widely used as a biomarker in crucial medical applications as well as for human age prediction of very high accuracy. This biomarker is based on the methylation status of several hundred CpG sites. In a recent line of publications we have adapted a versatile concept from evolutionary biology - the Universal Pacemaker (UPM) - to the setting of epigenetic aging and denoted it the Epigenetic PaceMaker (EPM). The EPM, as opposed to other epigenetic clocks, is not confined to specific pattern of aging, and the epigenetic age of the individual is inferred independently of other individuals. This allows an explicit modeling of aging trends, in particular non linear relationship between chronological and epigenetic age. In one of these recent works, we have presented an algorithmic improvement based on a two-step conditional expectation maximization (CEM) algorithm to arrive at a critical point on the likelihood surface. The algorithm alternates between a time step and a site step while advancing on the likelihood surface.ResultsHere we introduce non trivial improvements to these steps that are essential for analyzing data sets of realistic magnitude in a manageable time and space. These structural improvements are based on insights from linear algebra and symbolic algebra tools, providing us greater understanding of the degeneracy of the complex problem space. This understanding in turn, leads to the complete elimination of the bottleneck of cumbersome matrix multiplication and inversion, yielding a fast closed form solution in both steps of the CEM.In the experimental results part, we compare the CEM algorithm over several data sets and demonstrate the speedup obtained by the closed form solutions. Our results support the theoretical analysis of this improvement.ConclusionsThese improvements enable us to increase substantially the scale of inputs analyzed by the method, allowing us to apply the new approach to data sets that could not be analyzed before.

Automated localization of whales in coastal fjords.

Localization and tracking of vocalizing marine mammals are powerful tools for understanding and mitigating the impacts of anthropogenic stressors such as vessel noise on habitat use of cetaceans. A large-aperture hydrophone network has been installed in the Kitimat Fjord System, an ecologically, culturally, and economically valued marine environment in northern British Columbia, Canada. This network consists of four synchronized bottom-mounted hydrophones that permanently record and radio-transmit data to a land-based laboratory. An automated system has been developed which includes routines to localize transient bio-acoustic signals from three or more streaming hydrophones in near real-time. These routines comprise the correlation of hydrophone signals, the construction of a time lag model, and signal localization and error estimation from a spatial likelihood surface. The localization method was tested experimentally and subsequently applied to vocalizations from humpback whales, fin whales, and killer whales. Refractive and reflective sound propagation effects in the confined fjords are assessed using ray tracing propagation models. Automated localization results are compared to ground-truth data and shown to provide good accuracy.

Fast and Accurate Binary Response Mixed Model Analysis via Expectation Propagation

Expectation propagation is a general prescription for approximation of integrals in statistical inference problems. Its literature is mainly concerned with Bayesian inference scenarios. However, expectation propagation can also be used to approximate integrals arising in frequentist statistical inference. We focus on likelihood-based inference for binary response mixed models and show that fast and accurate quadrature-free inference can be realized for the probit link case with multivariate random effects and higher levels of nesting. The approach is supported by asymptotic calculations in which expectation propagation is seen to provide consistent estimation of the exact likelihood surface. Numerical studies reveal the availability of fast, highly accurate and scalable methodology for binary mixed model analysis.

Coalescence in the diffusion limit of a Bienaymé–Galton–Watson branching process

We consider the problem of estimating the elapsed time since the most recent common ancestor of a finite random sample drawn from a population which has evolved through a Bienaymé–Galton–Watson branching process. More specifically, we are interested in the diffusion limit appropriate to a supercritical process in the near-critical limit evolving over a large number of time steps. Our approach differs from earlier analyses in that we assume the only known information is the mean and variance of the number of offspring per parent, the observed total population size at the time of sampling, and the size of the sample. We obtain a formula for the probability that a finite random sample of the population is descended from a single ancestor in the initial population, and derive a confidence interval for the initial population size in terms of the final population size and the time since initiating the process. We also determine a joint likelihood surface from which confidence regions can be determined for simultaneously estimating two parameters, (1) the population size at the time of the most recent common ancestor, and (2) the time elapsed since the existence of the most recent common ancestor.